stanEx2
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}regression
normal linear model
ex5-1.stan
normal regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~normal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=normal_rng(m1[i],s);
}
}n=30
mdl=cmdstan_model('./ex5-1.stan')single regression
tb=tibble(x=runif(n,0,9),
y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)
qplot(data=tb,x,y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -8002.76
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -12.9393 0.000249625 0.000698878 1 1 52
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.mle## variable estimate
## lp__ -12.94
## b[1] 0.05
## b[2] 0.99
## s 0.93
## y1[1] 2.68
## y1[2] 3.09
## y1[3] 3.94
## y1[4] 5.56
## y1[5] 5.71
## y1[6] 5.54
## y1[7] 4.13
## y1[8] 3.11
## y1[9] 0.37
## y1[10] 8.16
## y1[11] 6.22
## y1[12] 4.65
## y1[13] 0.93
## y1[14] 2.82
## y1[15] 2.32
## y1[16] 3.06
## y1[17] 6.35
## y1[18] 9.67
## y1[19] 9.02
## y1[20] 3.89
## y1[21] 6.47
## y1[22] 5.50
## y1[23] 0.79
## y1[24] 2.50
## y1[25] 5.07
## y1[26] 0.58
## y1[27] 9.13
## y1[28] 0.87
## y1[29] 2.85
## y1[30] 2.76
## m1[1] 0.68
## m1[2] 2.84
## m1[3] 4.15
## m1[4] 5.05
## m1[5] 6.02
## m1[6] 5.07
## m1[7] 3.79
## m1[8] 0.89
## m1[9] 1.51
## m1[10] 8.86
## m1[11] 5.30
## m1[12] 4.44
## m1[13] 0.95
## m1[14] 3.59
## m1[15] 2.09
## m1[16] 2.28
## m1[17] 7.40
## m1[18] 8.43
## m1[19] 8.64
## m1[20] 2.88
## m1[21] 6.00
## m1[22] 5.47
## m1[23] 0.96
## m1[24] 2.81
## m1[25] 5.82
## m1[26] 1.49
## m1[27] 8.31
## m1[28] 1.27
## m1[29] 1.41
## m1[30] 4.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.62 -14.28 1.27 1.06 -17.12 -13.20 1.01 849 1150
## b[1] 0.07 0.08 0.36 0.36 -0.53 0.65 1.01 396 545
## b[2] 0.99 0.99 0.07 0.08 0.87 1.11 1.01 459 778
## s 1.02 1.00 0.15 0.14 0.81 1.28 1.00 1039 948
## y1[1] 0.68 0.65 1.04 1.03 -0.92 2.39 1.00 1794 1922
## y1[2] 2.86 2.86 1.04 1.04 1.16 4.53 1.00 1981 1686
## y1[3] 4.17 4.18 1.04 1.04 2.49 5.84 1.00 1965 1794
## y1[4] 5.06 5.08 1.02 1.00 3.41 6.78 1.00 1947 1933
## y1[5] 6.06 6.07 1.06 1.05 4.37 7.76 1.00 1780 1755
## y1[6] 5.10 5.11 1.06 1.03 3.36 6.85 1.00 1979 1855
## y1[7] 3.79 3.79 1.07 1.01 2.08 5.63 1.00 1831 1931
## y1[8] 0.92 0.94 1.06 1.03 -0.88 2.69 1.00 1515 1848
## y1[9] 1.51 1.49 1.03 1.01 -0.19 3.22 1.00 1686 1593
## y1[10] 8.90 8.88 1.12 1.06 7.08 10.71 1.00 2112 1759
## y1[11] 5.31 5.32 1.03 1.02 3.62 7.02 1.00 1789 1864
## y1[12] 4.45 4.47 1.02 0.94 2.79 6.16 1.00 1832 1601
## y1[13] 0.96 0.96 1.09 1.06 -0.82 2.73 1.00 1639 1775
## y1[14] 3.62 3.62 1.05 1.05 1.90 5.35 1.00 1919 1846
## y1[15] 2.06 2.05 1.05 0.97 0.28 3.80 1.00 1746 1879
## y1[16] 2.32 2.33 1.05 1.01 0.60 4.05 1.00 2088 1981
## y1[17] 7.41 7.38 1.08 1.06 5.64 9.18 1.00 1989 1971
## y1[18] 8.43 8.40 1.09 1.07 6.72 10.18 1.00 2102 1917
## y1[19] 8.65 8.64 1.11 1.09 6.91 10.51 1.00 1850 1863
## y1[20] 2.87 2.90 1.04 0.99 1.15 4.55 1.00 1928 1890
## y1[21] 5.97 5.96 1.06 1.01 4.21 7.72 1.00 1954 1859
## y1[22] 5.54 5.55 1.05 0.98 3.81 7.23 1.00 2184 1957
## y1[23] 0.96 0.99 1.08 1.05 -0.83 2.72 1.00 1808 1717
## y1[24] 2.81 2.83 1.07 0.99 1.07 4.61 1.00 1748 1984
## y1[25] 5.82 5.82 1.06 1.02 4.13 7.54 1.00 1925 1742
## y1[26] 1.50 1.49 1.05 1.04 -0.25 3.21 1.00 1397 1597
## y1[27] 8.31 8.29 1.09 1.07 6.49 10.08 1.00 1848 2104
## y1[28] 1.25 1.25 1.08 1.05 -0.58 3.02 1.00 1815 1931
## y1[29] 1.42 1.44 1.07 1.04 -0.39 3.18 1.00 1660 1913
## y1[30] 4.00 4.01 1.08 1.08 2.19 5.74 1.00 1836 1903
## m1[1] 0.70 0.71 0.32 0.32 0.16 1.21 1.01 410 474
## m1[2] 2.85 2.85 0.21 0.21 2.50 3.20 1.01 632 983
## m1[3] 4.16 4.15 0.19 0.19 3.85 4.47 1.00 1139 1265
## m1[4] 5.05 5.05 0.20 0.20 4.73 5.39 1.00 1730 1376
## m1[5] 6.02 6.01 0.24 0.23 5.65 6.40 1.00 1711 1540
## m1[6] 5.08 5.07 0.20 0.20 4.75 5.41 1.00 1738 1376
## m1[7] 3.80 3.80 0.19 0.19 3.49 4.12 1.00 940 1286
## m1[8] 0.91 0.91 0.31 0.31 0.39 1.39 1.01 417 477
## m1[9] 1.53 1.54 0.27 0.27 1.08 1.97 1.01 449 540
## m1[10] 8.86 8.86 0.40 0.40 8.20 9.51 1.01 798 1184
## m1[11] 5.30 5.29 0.21 0.20 4.96 5.65 1.00 1794 1400
## m1[12] 4.44 4.44 0.19 0.19 4.13 4.76 1.00 1335 1302
## m1[13] 0.96 0.97 0.30 0.30 0.45 1.45 1.01 420 477
## m1[14] 3.60 3.59 0.19 0.19 3.28 3.92 1.00 855 1249
## m1[15] 2.10 2.11 0.24 0.24 1.70 2.50 1.01 498 639
## m1[16] 2.29 2.29 0.23 0.24 1.90 2.67 1.01 522 726
## m1[17] 7.40 7.40 0.31 0.31 6.90 7.90 1.00 1032 1393
## m1[18] 8.43 8.43 0.37 0.37 7.82 9.03 1.01 847 1180
## m1[19] 8.64 8.64 0.39 0.39 8.00 9.26 1.01 822 1166
## m1[20] 2.89 2.89 0.21 0.21 2.54 3.24 1.01 640 982
## m1[21] 6.00 6.00 0.24 0.23 5.63 6.38 1.00 1717 1540
## m1[22] 5.48 5.47 0.21 0.21 5.14 5.83 1.00 1823 1421
## m1[23] 0.97 0.98 0.30 0.30 0.47 1.45 1.01 421 455
## m1[24] 2.82 2.81 0.21 0.21 2.46 3.17 1.01 622 907
## m1[25] 5.82 5.81 0.23 0.22 5.46 6.20 1.00 1778 1502
## m1[26] 1.50 1.51 0.27 0.27 1.05 1.94 1.01 448 540
## m1[27] 8.31 8.30 0.36 0.36 7.71 8.90 1.01 865 1154
## m1[28] 1.29 1.30 0.28 0.28 0.82 1.75 1.01 434 539
## m1[29] 1.42 1.43 0.28 0.27 0.96 1.87 1.01 441 507
## m1[30] 4.01 4.00 0.19 0.19 3.70 4.32 1.00 1050 1266
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
multiple regression
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -77.1658
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -13.335 0.000299111 0.00224091 1 1 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.34
## b[1] 0.27
## b[2] 1.04
## b[3] -1.06
## s 0.95
## y1[1] 4.97
## y1[2] -2.17
## y1[3] -0.15
## y1[4] 0.72
## y1[5] -6.24
## y1[6] 0.29
## y1[7] 2.00
## y1[8] -2.90
## y1[9] -6.79
## y1[10] -0.56
## y1[11] -1.64
## y1[12] 1.28
## y1[13] 5.10
## y1[14] 2.53
## y1[15] 0.93
## y1[16] 0.31
## y1[17] 1.55
## y1[18] 4.51
## y1[19] -3.39
## y1[20] 1.61
## y1[21] 3.69
## y1[22] 5.02
## y1[23] -1.66
## y1[24] 1.26
## y1[25] -1.39
## y1[26] -5.03
## y1[27] -1.85
## y1[28] 5.50
## y1[29] 7.91
## y1[30] 2.80
## m1[1] 5.18
## m1[2] -2.42
## m1[3] -0.93
## m1[4] -0.14
## m1[5] -6.83
## m1[6] 1.33
## m1[7] 2.48
## m1[8] -2.71
## m1[9] -6.60
## m1[10] -1.29
## m1[11] -0.59
## m1[12] 1.12
## m1[13] 3.26
## m1[14] 1.60
## m1[15] 0.39
## m1[16] 1.86
## m1[17] 0.32
## m1[18] 4.74
## m1[19] -3.54
## m1[20] -1.26
## m1[21] 3.51
## m1[22] 4.28
## m1[23] -2.33
## m1[24] 1.94
## m1[25] -1.03
## m1[26] -4.36
## m1[27] -2.96
## m1[28] 6.75
## m1[29] 7.23
## m1[30] 3.18mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.58 -15.19 1.58 1.31 -18.60 -13.83 1.01 635 906
## b[1] 0.28 0.29 0.41 0.42 -0.40 0.95 1.01 403 426
## b[2] 1.04 1.04 0.07 0.07 0.92 1.16 1.00 1226 1018
## b[3] -1.06 -1.06 0.07 0.07 -1.18 -0.93 1.00 740 1105
## s 1.05 1.03 0.16 0.15 0.82 1.34 1.00 941 1167
## y1[1] 5.16 5.14 1.09 1.03 3.41 6.90 1.00 2029 1710
## y1[2] -2.40 -2.43 1.14 1.07 -4.24 -0.49 1.00 1895 1802
## y1[3] -0.94 -0.92 1.12 1.06 -2.75 0.92 1.00 1413 1880
## y1[4] -0.12 -0.13 1.09 1.09 -1.85 1.66 1.00 1711 1951
## y1[5] -6.81 -6.83 1.15 1.13 -8.59 -4.86 1.00 2099 1540
## y1[6] 1.33 1.33 1.11 1.07 -0.53 3.13 1.00 1954 1548
## y1[7] 2.47 2.47 1.13 1.09 0.63 4.32 1.00 1628 1900
## y1[8] -2.72 -2.73 1.10 1.08 -4.43 -0.95 1.00 2006 1785
## y1[9] -6.56 -6.55 1.12 1.09 -8.45 -4.72 1.00 2045 1934
## y1[10] -1.26 -1.23 1.08 1.06 -3.01 0.48 1.00 1782 1838
## y1[11] -0.57 -0.57 1.08 1.09 -2.38 1.18 1.00 1926 1737
## y1[12] 1.07 1.06 1.13 1.11 -0.72 2.94 1.00 1541 1843
## y1[13] 3.24 3.25 1.09 1.07 1.44 5.03 1.00 1910 1805
## y1[14] 1.55 1.55 1.08 1.07 -0.19 3.36 1.00 1950 1968
## y1[15] 0.38 0.37 1.08 1.05 -1.41 2.17 1.00 2062 1961
## y1[16] 1.89 1.91 1.11 1.05 0.07 3.68 1.00 1537 1767
## y1[17] 0.35 0.35 1.16 1.13 -1.52 2.24 1.00 1697 1786
## y1[18] 4.71 4.74 1.10 1.06 2.93 6.46 1.00 1778 1778
## y1[19] -3.50 -3.52 1.10 1.12 -5.29 -1.70 1.00 1873 1840
## y1[20] -1.27 -1.30 1.12 1.10 -3.10 0.61 1.00 2058 2020
## y1[21] 3.55 3.56 1.12 1.11 1.69 5.36 1.00 1792 1895
## y1[22] 4.25 4.28 1.12 1.05 2.34 6.09 1.00 1865 1822
## y1[23] -2.33 -2.34 1.10 1.08 -4.12 -0.53 1.00 2027 1930
## y1[24] 1.92 1.95 1.11 1.07 0.15 3.74 1.00 1736 1774
## y1[25] -1.08 -1.09 1.10 1.08 -2.87 0.73 1.00 1790 1875
## y1[26] -4.35 -4.35 1.13 1.09 -6.22 -2.54 1.00 1797 1857
## y1[27] -2.98 -2.97 1.09 1.07 -4.79 -1.24 1.00 2116 1916
## y1[28] 6.72 6.70 1.12 1.08 4.97 8.59 1.00 2012 1881
## y1[29] 7.23 7.19 1.15 1.12 5.42 9.17 1.00 1761 1684
## y1[30] 3.17 3.18 1.10 1.05 1.41 5.04 1.00 1867 1871
## m1[1] 5.17 5.17 0.34 0.34 4.62 5.72 1.00 1674 1279
## m1[2] -2.42 -2.42 0.37 0.35 -3.00 -1.81 1.01 1079 1195
## m1[3] -0.93 -0.91 0.35 0.35 -1.50 -0.37 1.01 441 451
## m1[4] -0.14 -0.13 0.29 0.29 -0.62 0.32 1.01 437 467
## m1[5] -6.83 -6.84 0.46 0.45 -7.56 -6.07 1.01 1993 1324
## m1[6] 1.33 1.33 0.33 0.33 0.79 1.84 1.01 1256 1088
## m1[7] 2.48 2.49 0.34 0.35 1.92 3.03 1.01 428 554
## m1[8] -2.70 -2.71 0.26 0.25 -3.12 -2.27 1.00 1478 1354
## m1[9] -6.59 -6.60 0.45 0.43 -7.31 -5.87 1.01 1997 1334
## m1[10] -1.28 -1.27 0.24 0.24 -1.70 -0.90 1.01 614 721
## m1[11] -0.59 -0.58 0.22 0.21 -0.95 -0.23 1.01 1970 1414
## m1[12] 1.11 1.11 0.39 0.38 0.46 1.71 1.01 1007 1045
## m1[13] 3.25 3.27 0.32 0.31 2.71 3.76 1.01 1735 1326
## m1[14] 1.60 1.60 0.22 0.23 1.21 1.94 1.01 1723 1468
## m1[15] 0.38 0.39 0.26 0.26 -0.04 0.78 1.01 1680 1298
## m1[16] 1.86 1.87 0.36 0.37 1.26 2.45 1.01 406 459
## m1[17] 0.32 0.32 0.43 0.41 -0.40 0.99 1.01 871 937
## m1[18] 4.73 4.73 0.34 0.33 4.19 5.27 1.00 1924 1518
## m1[19] -3.53 -3.54 0.31 0.30 -4.05 -3.02 1.00 995 1248
## m1[20] -1.26 -1.26 0.22 0.21 -1.62 -0.90 1.01 1676 1222
## m1[21] 3.50 3.51 0.36 0.35 2.90 4.06 1.01 1511 1267
## m1[22] 4.28 4.28 0.32 0.31 3.74 4.81 1.00 652 967
## m1[23] -2.32 -2.32 0.26 0.25 -2.76 -1.89 1.01 2143 1391
## m1[24] 1.94 1.94 0.35 0.35 1.36 2.51 1.01 410 490
## m1[25] -1.03 -1.01 0.36 0.36 -1.61 -0.45 1.01 442 425
## m1[26] -4.35 -4.35 0.36 0.36 -4.94 -3.76 1.00 914 1268
## m1[27] -2.96 -2.96 0.27 0.26 -3.40 -2.50 1.00 1998 1402
## m1[28] 6.74 6.74 0.41 0.40 6.09 7.41 1.00 1605 1304
## m1[29] 7.22 7.22 0.44 0.43 6.53 7.94 1.00 1098 1248
## m1[30] 3.18 3.18 0.27 0.27 2.72 3.62 1.00 640 949
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANOVA
tb=tibble(c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)
qplot(data=tb,c,y,geom='boxplot')
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -41.02
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -15.6149 0.000624738 0.000527558 1 1 20
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -15.61
## b[1] -0.23
## b[2] 2.36
## b[3] -1.82
## s 1.02
## y1[1] 2.84
## y1[2] -1.61
## y1[3] -1.08
## y1[4] 1.85
## y1[5] -0.88
## y1[6] 3.70
## y1[7] -0.40
## y1[8] -1.28
## y1[9] 0.40
## y1[10] 3.53
## y1[11] -0.22
## y1[12] -2.51
## y1[13] -0.41
## y1[14] 1.97
## y1[15] 2.60
## y1[16] 2.48
## y1[17] -0.94
## y1[18] 3.08
## y1[19] -2.31
## y1[20] 0.47
## y1[21] -2.97
## y1[22] -2.25
## y1[23] -2.68
## y1[24] -1.02
## y1[25] 1.34
## y1[26] -0.44
## y1[27] 0.33
## y1[28] 0.18
## y1[29] -1.39
## y1[30] 2.63
## m1[1] 2.13
## m1[2] -2.05
## m1[3] -0.23
## m1[4] 2.13
## m1[5] -0.23
## m1[6] 2.13
## m1[7] -0.23
## m1[8] -2.05
## m1[9] -0.23
## m1[10] 2.13
## m1[11] -0.23
## m1[12] -2.05
## m1[13] -0.23
## m1[14] 2.13
## m1[15] 2.13
## m1[16] 2.13
## m1[17] -0.23
## m1[18] 2.13
## m1[19] -2.05
## m1[20] -0.23
## m1[21] -2.05
## m1[22] -2.05
## m1[23] -2.05
## m1[24] -2.05
## m1[25] 2.13
## m1[26] -0.23
## m1[27] -0.23
## m1[28] -0.23
## m1[29] -2.05
## m1[30] 2.13mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.78 -17.43 1.55 1.31 -20.81 -15.97 1.00 742 1057
## b[1] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## b[2] 2.35 2.34 0.49 0.48 1.57 3.15 1.00 718 781
## b[3] -1.83 -1.83 0.55 0.52 -2.75 -0.92 1.00 635 740
## s 1.13 1.11 0.17 0.16 0.90 1.42 1.00 1894 1741
## y1[1] 2.11 2.14 1.21 1.17 0.15 4.07 1.00 1957 1893
## y1[2] -2.09 -2.07 1.20 1.16 -4.05 -0.17 1.00 1709 1833
## y1[3] -0.20 -0.18 1.21 1.18 -2.25 1.77 1.00 1740 1921
## y1[4] 2.15 2.17 1.22 1.20 0.14 4.15 1.00 1675 1536
## y1[5] -0.19 -0.20 1.20 1.19 -2.11 1.81 1.00 1596 1697
## y1[6] 2.12 2.11 1.21 1.22 0.14 4.17 1.00 1838 1999
## y1[7] -0.22 -0.18 1.17 1.16 -2.19 1.67 1.00 1862 1735
## y1[8] -2.07 -2.11 1.22 1.17 -4.07 -0.05 1.00 1883 1860
## y1[9] -0.21 -0.20 1.19 1.21 -2.19 1.76 1.00 1426 1901
## y1[10] 2.16 2.16 1.23 1.20 0.16 4.11 1.00 1743 1811
## y1[11] -0.24 -0.24 1.20 1.12 -2.15 1.76 1.00 1508 1956
## y1[12] -2.06 -2.04 1.21 1.19 -4.08 -0.10 1.00 1892 1700
## y1[13] -0.25 -0.25 1.18 1.12 -2.17 1.75 1.00 1747 1973
## y1[14] 2.11 2.13 1.18 1.17 0.15 4.01 1.00 2076 1973
## y1[15] 2.12 2.10 1.17 1.13 0.20 4.02 1.00 2012 1994
## y1[16] 2.15 2.15 1.18 1.14 0.22 4.13 1.00 1913 1643
## y1[17] -0.23 -0.22 1.19 1.18 -2.22 1.71 1.00 1751 1915
## y1[18] 2.14 2.11 1.24 1.26 0.12 4.19 1.00 1924 1855
## y1[19] -2.05 -2.04 1.19 1.16 -4.01 -0.11 1.00 1835 1806
## y1[20] -0.21 -0.21 1.19 1.17 -2.12 1.79 1.00 1812 1774
## y1[21] -2.06 -2.05 1.24 1.18 -4.07 0.00 1.00 1849 1744
## y1[22] -2.08 -2.10 1.21 1.13 -4.07 -0.09 1.00 1695 1938
## y1[23] -2.08 -2.06 1.22 1.18 -4.08 -0.07 1.00 1676 1680
## y1[24] -2.05 -2.03 1.20 1.19 -4.05 -0.05 1.00 1855 2013
## y1[25] 2.12 2.11 1.21 1.23 0.15 4.10 1.00 1975 1931
## y1[26] -0.22 -0.18 1.19 1.11 -2.26 1.70 1.00 1784 1899
## y1[27] -0.24 -0.21 1.19 1.18 -2.23 1.68 1.00 1864 1846
## y1[28] -0.25 -0.24 1.21 1.18 -2.22 1.77 1.00 1737 1910
## y1[29] -2.02 -2.02 1.20 1.19 -3.96 -0.03 1.00 1832 1971
## y1[30] 2.12 2.08 1.21 1.15 0.15 4.22 1.00 1867 1963
## m1[1] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[2] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[3] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[4] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[5] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[6] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[7] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[8] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[9] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[10] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[11] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[12] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[13] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[14] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[15] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[16] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[17] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[18] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[19] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[20] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[21] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[22] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[23] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[24] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[25] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
## m1[26] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[27] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[28] -0.22 -0.21 0.35 0.33 -0.79 0.36 1.00 709 916
## m1[29] -2.05 -2.05 0.40 0.38 -2.71 -1.41 1.00 1205 1120
## m1[30] 2.13 2.13 0.36 0.35 1.53 2.70 1.00 1404 1387
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANCOVA
tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))
f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)
lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])
qplot(data=tb,x,y,shape=c,size=I(2))
plot(tb$x,tb$y,col=1+lv[factor(tb$c)])
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -5073.62
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 38 -14.9114 4.98066e-05 0.00237982 0.2792 1 69
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.91
## b[1] -0.25
## b[2] 1.03
## b[3] 2.00
## b[4] -2.81
## s 1.00
## y1[1] 8.11
## y1[2] -0.56
## y1[3] 1.60
## y1[4] 4.14
## y1[5] 8.45
## y1[6] 6.87
## y1[7] 4.97
## y1[8] 2.45
## y1[9] 1.75
## y1[10] 8.35
## y1[11] 4.08
## y1[12] 4.12
## y1[13] 5.94
## y1[14] -2.97
## y1[15] 6.52
## y1[16] 1.97
## y1[17] 10.81
## y1[18] 6.85
## y1[19] 1.38
## y1[20] 2.96
## y1[21] 4.67
## y1[22] 5.19
## y1[23] 12.04
## y1[24] 4.15
## y1[25] 1.46
## y1[26] 5.92
## y1[27] 2.57
## y1[28] 3.58
## y1[29] 1.33
## y1[30] 3.03
## m1[1] 9.23
## m1[2] -2.68
## m1[3] 1.67
## m1[4] 4.12
## m1[5] 8.70
## m1[6] 8.27
## m1[7] 5.42
## m1[8] 1.85
## m1[9] 1.85
## m1[10] 9.26
## m1[11] 2.42
## m1[12] 4.36
## m1[13] 7.03
## m1[14] -2.14
## m1[15] 6.39
## m1[16] 1.65
## m1[17] 8.92
## m1[18] 7.32
## m1[19] 1.90
## m1[20] 1.51
## m1[21] 6.72
## m1[22] 5.22
## m1[23] 10.57
## m1[24] 3.17
## m1[25] 0.32
## m1[26] 5.47
## m1[27] 2.65
## m1[28] 5.07
## m1[29] 1.26
## m1[30] 2.92mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.69 -17.31 1.76 1.55 -20.89 -15.55 1.00 609 1054
## b[1] -0.25 -0.27 0.46 0.43 -0.98 0.53 1.00 726 988
## b[2] 1.03 1.03 0.08 0.08 0.90 1.16 1.00 1598 1659
## b[3] 2.00 2.02 0.49 0.47 1.18 2.82 1.00 872 906
## b[4] -2.83 -2.82 0.60 0.57 -3.82 -1.81 1.01 471 821
## s 1.12 1.11 0.17 0.16 0.89 1.42 1.00 1720 1177
## y1[1] 9.25 9.25 1.18 1.19 7.31 11.18 1.00 2239 1904
## y1[2] -2.64 -2.63 1.27 1.23 -4.69 -0.57 1.00 1880 1930
## y1[3] 1.67 1.66 1.19 1.19 -0.25 3.61 1.00 1726 2013
## y1[4] 4.12 4.10 1.20 1.19 2.17 6.07 1.00 2015 1931
## y1[5] 8.69 8.71 1.20 1.12 6.65 10.62 1.00 1946 1918
## y1[6] 8.25 8.26 1.17 1.09 6.32 10.13 1.00 1852 2014
## y1[7] 5.43 5.42 1.26 1.19 3.38 7.46 1.00 1598 1878
## y1[8] 1.86 1.84 1.26 1.24 -0.19 3.97 1.00 1918 1979
## y1[9] 1.81 1.82 1.19 1.10 -0.12 3.75 1.00 1832 1631
## y1[10] 9.24 9.24 1.19 1.21 7.31 11.15 1.00 1994 1881
## y1[11] 2.45 2.46 1.25 1.20 0.37 4.53 1.00 1876 1854
## y1[12] 4.30 4.30 1.17 1.13 2.37 6.24 1.00 1795 2015
## y1[13] 7.03 7.02 1.16 1.13 5.16 8.97 1.00 1875 1974
## y1[14] -2.18 -2.19 1.22 1.22 -4.18 -0.16 1.00 1796 2031
## y1[15] 6.38 6.39 1.22 1.16 4.31 8.39 1.00 1818 1685
## y1[16] 1.68 1.70 1.22 1.18 -0.35 3.65 1.00 1813 1937
## y1[17] 8.93 8.94 1.19 1.13 7.03 10.90 1.00 1608 1743
## y1[18] 7.30 7.30 1.14 1.08 5.42 9.24 1.00 1860 1927
## y1[19] 1.90 1.87 1.21 1.19 0.01 3.87 1.00 1957 1684
## y1[20] 1.49 1.51 1.22 1.21 -0.52 3.43 1.00 1685 1932
## y1[21] 6.74 6.73 1.24 1.17 4.72 8.85 1.00 1863 1933
## y1[22] 5.21 5.22 1.22 1.17 3.14 7.21 1.00 1712 1540
## y1[23] 10.59 10.61 1.23 1.15 8.52 12.57 1.00 1669 1582
## y1[24] 3.18 3.19 1.23 1.18 1.06 5.15 1.00 1960 1932
## y1[25] 0.31 0.32 1.25 1.20 -1.72 2.33 1.00 1531 1836
## y1[26] 5.49 5.49 1.16 1.17 3.63 7.43 1.00 2064 2000
## y1[27] 2.65 2.64 1.16 1.11 0.76 4.58 1.00 1795 1951
## y1[28] 5.03 5.04 1.18 1.14 3.07 6.98 1.00 2029 1896
## y1[29] 1.24 1.27 1.18 1.18 -0.72 3.22 1.00 1753 1912
## y1[30] 2.93 2.94 1.21 1.17 0.94 4.87 1.00 2043 1835
## m1[1] 9.24 9.24 0.38 0.37 8.62 9.85 1.00 1620 1447
## m1[2] -2.70 -2.70 0.53 0.50 -3.54 -1.82 1.00 928 1062
## m1[3] 1.67 1.64 0.39 0.37 1.04 2.30 1.00 643 877
## m1[4] 4.12 4.12 0.38 0.38 3.49 4.74 1.00 1435 1079
## m1[5] 8.70 8.70 0.36 0.35 8.11 9.27 1.00 1597 1470
## m1[6] 8.27 8.27 0.34 0.33 7.71 8.83 1.00 1578 1469
## m1[7] 5.41 5.42 0.61 0.58 4.42 6.39 1.00 1051 1465
## m1[8] 1.84 1.85 0.50 0.50 1.04 2.65 1.00 1451 881
## m1[9] 1.85 1.83 0.38 0.36 1.23 2.48 1.00 632 866
## m1[10] 9.26 9.27 0.38 0.37 8.64 9.88 1.00 1622 1447
## m1[11] 2.40 2.41 0.49 0.49 1.61 3.18 1.01 883 1215
## m1[12] 4.36 4.35 0.38 0.37 3.74 4.99 1.00 705 851
## m1[13] 7.03 7.03 0.32 0.31 6.52 7.56 1.00 1511 1567
## m1[14] -2.15 -2.15 0.51 0.50 -2.97 -1.31 1.00 897 988
## m1[15] 6.40 6.40 0.44 0.42 5.65 7.11 1.00 1050 1270
## m1[16] 1.64 1.65 0.48 0.47 0.86 2.41 1.01 858 1251
## m1[17] 8.92 8.93 0.36 0.36 8.32 9.51 1.00 1607 1404
## m1[18] 7.32 7.32 0.32 0.31 6.79 7.84 1.00 1527 1541
## m1[19] 1.90 1.88 0.38 0.36 1.28 2.52 1.00 631 806
## m1[20] 1.50 1.51 0.48 0.47 0.72 2.26 1.01 854 1348
## m1[21] 6.72 6.73 0.46 0.44 5.95 7.45 1.00 1090 1263
## m1[22] 5.22 5.22 0.40 0.38 4.57 5.89 1.00 786 1087
## m1[23] 10.58 10.58 0.44 0.44 9.84 11.29 1.00 1655 1516
## m1[24] 3.17 3.17 0.43 0.43 2.47 3.86 1.00 1434 961
## m1[25] 0.31 0.29 0.43 0.41 -0.38 1.05 1.00 705 958
## m1[26] 5.47 5.47 0.33 0.33 4.90 6.01 1.00 1455 1212
## m1[27] 2.65 2.63 0.37 0.36 2.05 3.27 1.00 626 882
## m1[28] 5.07 5.07 0.35 0.34 4.49 5.63 1.00 1445 1195
## m1[29] 1.26 1.24 0.40 0.38 0.62 1.93 1.00 666 889
## m1[30] 2.92 2.93 0.44 0.44 2.20 3.63 1.00 1436 966
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
interaction of variable
n=50
mdl=cmdstan_model('./ex5-1.stan')
tb=tibble(x=runif(n,-3,3),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,x,y,col=ca),
qplot(data=tb,x,y,col=cb),ncol=2)
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)
fn(f0)## Initial log joint probability = -2013.08
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -20.4747 7.45039e-05 0.000373311 0.8734 0.8734 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.59 -22.30 1.43 1.26 -25.35 -20.94 1.00 811 1095
## b[1] 0.33 0.33 0.17 0.16 0.04 0.61 1.00 765 972
## b[2] 1.22 1.22 0.08 0.08 1.09 1.36 1.00 2614 1715
## b[3] 0.44 0.44 0.29 0.27 -0.03 0.93 1.00 630 843
## s 0.96 0.96 0.10 0.10 0.81 1.15 1.00 1850 1570
## y1[1] 1.31 1.27 0.97 0.97 -0.29 2.90 1.00 2035 2057
## y1[2] 1.07 1.04 0.97 0.94 -0.49 2.68 1.00 2021 1967
## y1[3] 3.16 3.16 0.99 0.95 1.54 4.75 1.00 1971 1854
## y1[4] 2.68 2.72 1.01 0.99 0.99 4.38 1.00 2033 1988
## y1[5] -2.98 -2.98 1.02 1.01 -4.60 -1.34 1.00 1845 1930
## y1[6] -1.62 -1.59 1.00 0.96 -3.32 0.05 1.00 1737 1912
## y1[7] 0.58 0.59 0.96 0.93 -0.94 2.18 1.00 1863 1893
## y1[8] -1.34 -1.33 1.00 0.98 -3.05 0.30 1.00 1954 1848
## y1[9] -0.85 -0.86 0.99 0.98 -2.45 0.76 1.00 1889 1817
## y1[10] 2.35 2.37 0.99 0.97 0.72 4.00 1.00 2050 1973
## y1[11] 2.37 2.39 1.01 1.00 0.72 4.04 1.00 2154 1907
## y1[12] 0.02 0.01 0.99 0.99 -1.60 1.65 1.00 1931 1971
## y1[13] -1.85 -1.87 1.02 1.04 -3.55 -0.14 1.00 2077 2055
## y1[14] 3.21 3.19 1.01 0.97 1.61 4.87 1.00 2107 2048
## y1[15] -3.27 -3.26 1.01 0.99 -5.01 -1.57 1.00 1923 1660
## y1[16] 0.44 0.45 0.97 0.94 -1.12 2.01 1.00 1888 1959
## y1[17] 2.79 2.78 1.01 1.01 1.16 4.44 1.00 1876 1883
## y1[18] -1.56 -1.58 0.99 0.97 -3.15 0.07 1.00 1901 1880
## y1[19] 3.33 3.35 1.01 1.02 1.67 5.00 1.00 1876 1822
## y1[20] 3.65 3.65 1.00 0.99 2.00 5.29 1.00 1931 1844
## y1[21] 2.17 2.17 0.98 0.99 0.57 3.77 1.00 1865 1959
## y1[22] -1.83 -1.85 1.03 0.99 -3.55 -0.13 1.00 1811 1933
## y1[23] 1.58 1.60 0.98 0.95 -0.03 3.23 1.00 1898 2103
## y1[24] 0.41 0.42 1.00 0.95 -1.24 2.02 1.00 1883 1904
## y1[25] 1.38 1.39 0.95 0.97 -0.23 2.91 1.00 1757 1772
## y1[26] -0.20 -0.19 1.00 0.99 -1.88 1.47 1.00 1829 1977
## y1[27] 2.10 2.14 0.98 0.97 0.43 3.68 1.00 1807 1970
## y1[28] 0.96 0.93 0.97 0.99 -0.65 2.55 1.00 2021 2000
## y1[29] -2.17 -2.19 1.02 1.02 -3.82 -0.43 1.00 1709 1790
## y1[30] 2.56 2.58 0.99 0.94 0.89 4.18 1.00 1719 1839
## y1[31] 2.06 2.04 1.00 0.95 0.41 3.71 1.00 1821 1931
## y1[32] 4.29 4.30 1.01 1.00 2.58 5.93 1.00 1864 1974
## y1[33] 3.94 3.96 1.03 1.02 2.20 5.58 1.00 1816 1700
## y1[34] 2.58 2.58 1.01 0.98 0.89 4.23 1.00 1460 2001
## y1[35] 1.14 1.16 1.00 0.94 -0.52 2.75 1.00 1970 1558
## y1[36] -0.20 -0.18 1.00 0.99 -1.88 1.45 1.00 1934 2010
## y1[37] -0.61 -0.63 0.96 0.92 -2.16 1.00 1.00 1975 1968
## y1[38] -1.51 -1.51 0.99 0.95 -3.14 0.11 1.00 1928 2054
## y1[39] 3.14 3.13 1.00 0.98 1.50 4.78 1.00 1807 1955
## y1[40] -2.97 -2.96 1.00 1.00 -4.56 -1.38 1.00 1854 1917
## y1[41] 2.88 2.93 1.01 1.00 1.20 4.54 1.00 1864 1819
## y1[42] 1.99 1.98 0.98 0.97 0.41 3.61 1.00 2006 1933
## y1[43] 0.37 0.39 1.00 0.97 -1.26 1.95 1.00 1894 1968
## y1[44] 3.92 3.93 1.02 1.02 2.17 5.54 1.00 2036 2011
## y1[45] -1.41 -1.42 1.01 0.99 -3.08 0.29 1.00 1837 1807
## y1[46] 2.57 2.58 1.00 0.99 0.99 4.24 1.00 1807 1852
## y1[47] -0.48 -0.51 1.00 0.98 -2.14 1.13 1.00 1797 2057
## y1[48] 1.64 1.64 0.98 0.92 -0.05 3.28 1.00 1663 1729
## y1[49] 2.88 2.88 1.01 0.98 1.23 4.51 1.00 1860 1783
## y1[50] -0.24 -0.25 0.99 0.94 -1.82 1.36 1.00 1818 1777
## m1[1] 1.27 1.28 0.17 0.17 0.99 1.57 1.00 817 1207
## m1[2] 1.07 1.08 0.17 0.16 0.79 1.36 1.00 793 1068
## m1[3] 3.17 3.17 0.26 0.25 2.73 3.59 1.00 1189 1457
## m1[4] 2.69 2.69 0.25 0.23 2.29 3.09 1.00 1109 1431
## m1[5] -2.94 -2.94 0.30 0.29 -3.43 -2.45 1.00 1374 1492
## m1[6] -1.63 -1.62 0.29 0.29 -2.10 -1.16 1.00 1301 1471
## m1[7] 0.62 0.62 0.17 0.16 0.34 0.90 1.00 762 1011
## m1[8] -1.36 -1.37 0.22 0.21 -1.74 -1.00 1.00 1039 989
## m1[9] -0.83 -0.83 0.20 0.19 -1.16 -0.50 1.00 930 1070
## m1[10] 2.35 2.35 0.20 0.20 2.02 2.70 1.00 1061 1330
## m1[11] 2.38 2.37 0.24 0.23 1.99 2.77 1.00 1062 1390
## m1[12] 0.02 0.02 0.18 0.17 -0.28 0.31 1.00 786 986
## m1[13] -1.84 -1.83 0.30 0.30 -2.33 -1.36 1.00 1340 1448
## m1[14] 3.21 3.21 0.24 0.23 2.82 3.60 1.00 1356 1595
## m1[15] -3.27 -3.27 0.32 0.31 -3.80 -2.76 1.00 1440 1547
## m1[16] 0.46 0.46 0.17 0.16 0.17 0.74 1.00 762 948
## m1[17] 2.79 2.79 0.22 0.21 2.44 3.16 1.00 1206 1460
## m1[18] -1.56 -1.57 0.23 0.22 -1.95 -1.19 1.00 1084 1097
## m1[19] 3.32 3.32 0.25 0.23 2.92 3.72 1.00 1396 1595
## m1[20] 3.66 3.66 0.26 0.25 3.24 4.09 1.00 1523 1631
## m1[21] 2.20 2.20 0.20 0.19 1.88 2.53 1.00 1016 1395
## m1[22] -1.83 -1.82 0.30 0.30 -2.32 -1.35 1.00 1339 1448
## m1[23] 1.57 1.58 0.18 0.17 1.28 1.88 1.00 859 1288
## m1[24] 0.41 0.41 0.17 0.16 0.12 0.69 1.00 762 998
## m1[25] 1.36 1.36 0.18 0.17 1.07 1.65 1.00 824 1271
## m1[26] -0.21 -0.21 0.24 0.24 -0.60 0.18 1.00 1063 1268
## m1[27] 2.10 2.11 0.19 0.19 1.79 2.43 1.00 988 1343
## m1[28] 0.97 0.96 0.22 0.22 0.61 1.34 1.00 975 1226
## m1[29] -2.21 -2.20 0.31 0.31 -2.73 -1.71 1.00 1412 1461
## m1[30] 2.59 2.59 0.21 0.20 2.25 2.95 1.00 1139 1440
## m1[31] 2.07 2.07 0.23 0.22 1.69 2.45 1.00 1029 1251
## m1[32] 4.31 4.32 0.31 0.30 3.79 4.81 1.00 1413 1554
## m1[33] 3.94 3.95 0.29 0.28 3.44 4.42 1.00 1337 1516
## m1[34] 2.57 2.57 0.21 0.20 2.23 2.93 1.00 1131 1369
## m1[35] 1.12 1.12 0.17 0.16 0.83 1.41 1.00 799 1115
## m1[36] -0.19 -0.19 0.18 0.17 -0.50 0.10 1.00 817 1025
## m1[37] -0.62 -0.62 0.19 0.18 -0.94 -0.30 1.00 893 1036
## m1[38] -1.51 -1.52 0.23 0.22 -1.90 -1.14 1.00 1073 1083
## m1[39] 3.19 3.19 0.26 0.25 2.75 3.62 1.00 1194 1505
## m1[40] -3.00 -3.00 0.30 0.29 -3.50 -2.51 1.00 1386 1447
## m1[41] 2.92 2.91 0.25 0.24 2.50 3.33 1.00 1146 1475
## m1[42] 1.98 1.98 0.23 0.22 1.61 2.36 1.00 1020 1220
## m1[43] 0.35 0.35 0.23 0.23 -0.02 0.73 1.00 1006 1233
## m1[44] 3.93 3.93 0.28 0.27 3.49 4.39 1.00 1626 1624
## m1[45] -1.40 -1.40 0.22 0.21 -1.77 -1.03 1.00 1046 989
## m1[46] 2.57 2.57 0.21 0.20 2.23 2.93 1.00 1131 1369
## m1[47] -0.50 -0.50 0.25 0.26 -0.91 -0.10 1.00 1101 1309
## m1[48] 1.62 1.62 0.22 0.22 1.26 1.99 1.00 991 1238
## m1[49] 2.89 2.89 0.23 0.21 2.53 3.26 1.00 1242 1482
## m1[50] -0.25 -0.25 0.18 0.17 -0.56 0.05 1.00 826 1003
fn(f1)## Initial log joint probability = -99.5729
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -14.4515 8.86349e-06 0.00140858 0.4152 0.4152 29
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.31 -16.93 1.78 1.45 -20.65 -15.26 1.00 776 694
## b[1] -0.04 -0.04 0.19 0.19 -0.36 0.27 1.00 1243 991
## b[2] 1.18 1.18 0.07 0.07 1.06 1.30 1.00 2515 1487
## b[3] 0.45 0.44 0.27 0.25 0.03 0.91 1.00 1002 709
## b[4] 0.87 0.87 0.24 0.24 0.47 1.27 1.00 1073 946
## s 0.87 0.86 0.10 0.09 0.73 1.04 1.00 2343 1624
## y1[1] 0.86 0.86 0.91 0.91 -0.58 2.31 1.00 1771 1817
## y1[2] 0.63 0.64 0.89 0.89 -0.85 2.14 1.00 2054 1972
## y1[3] 3.59 3.61 0.89 0.86 2.09 5.01 1.00 1662 1663
## y1[4] 2.22 2.19 0.91 0.90 0.78 3.76 1.00 1757 1954
## y1[5] -2.31 -2.31 0.93 0.89 -3.85 -0.78 1.00 2015 1830
## y1[6] -1.89 -1.90 0.93 0.90 -3.40 -0.40 1.00 1982 1877
## y1[7] 1.10 1.11 0.89 0.88 -0.32 2.59 1.00 1844 1811
## y1[8] -1.73 -1.75 0.88 0.83 -3.19 -0.24 1.00 1844 1877
## y1[9] -0.30 -0.30 0.92 0.90 -1.76 1.15 1.00 2100 1967
## y1[10] 1.87 1.86 0.88 0.87 0.41 3.29 1.00 1851 1859
## y1[11] 1.95 1.97 0.87 0.86 0.51 3.34 1.00 1977 1974
## y1[12] 0.48 0.49 0.91 0.86 -1.01 1.97 1.00 2018 2014
## y1[13] -2.11 -2.11 0.91 0.90 -3.55 -0.60 1.00 2120 1882
## y1[14] 2.72 2.73 0.91 0.88 1.16 4.20 1.00 1817 1895
## y1[15] -3.53 -3.54 0.91 0.91 -4.98 -2.04 1.00 1691 1808
## y1[16] 0.06 0.07 0.92 0.90 -1.51 1.61 1.00 1837 1912
## y1[17] 2.34 2.34 0.90 0.92 0.87 3.82 1.00 1603 1945
## y1[18] -1.00 -1.00 0.90 0.91 -2.43 0.47 1.00 1859 1943
## y1[19] 2.83 2.82 0.91 0.91 1.36 4.31 1.00 2012 1959
## y1[20] 4.04 4.01 0.91 0.90 2.57 5.57 1.00 1821 2046
## y1[21] 1.77 1.75 0.90 0.87 0.34 3.31 1.00 2099 1689
## y1[22] -2.12 -2.10 0.91 0.88 -3.65 -0.67 1.00 1849 1852
## y1[23] 2.04 2.03 0.90 0.87 0.51 3.46 1.00 1932 1939
## y1[24] 0.01 0.02 0.88 0.89 -1.41 1.46 1.00 1931 1918
## y1[25] 0.96 0.94 0.92 0.87 -0.55 2.44 1.00 1761 1820
## y1[26] -0.54 -0.54 0.92 0.89 -2.04 0.98 1.00 1797 1998
## y1[27] 1.68 1.66 0.91 0.89 0.20 3.21 1.00 2076 1796
## y1[28] 1.47 1.48 0.93 0.92 -0.06 3.00 1.00 1705 2012
## y1[29] -1.57 -1.55 0.93 0.92 -3.08 -0.06 1.00 2001 1955
## y1[30] 3.03 3.02 0.88 0.86 1.60 4.49 1.00 2062 1845
## y1[31] 2.51 2.51 0.92 0.92 1.00 4.03 1.00 1882 1969
## y1[32] 3.84 3.81 0.93 0.92 2.28 5.38 1.00 1895 1833
## y1[33] 4.35 4.37 0.91 0.91 2.86 5.81 1.00 1842 1933
## y1[34] 2.98 2.97 0.89 0.85 1.54 4.43 1.00 2005 1970
## y1[35] 1.59 1.57 0.89 0.89 0.19 3.07 1.00 2039 1932
## y1[36] 0.33 0.33 0.92 0.90 -1.23 1.81 1.00 1886 1699
## y1[37] -0.10 -0.09 0.89 0.87 -1.51 1.39 1.00 1840 1974
## y1[38] -1.84 -1.82 0.90 0.89 -3.37 -0.38 1.00 2117 1871
## y1[39] 3.65 3.66 0.91 0.91 2.16 5.17 1.00 1880 1915
## y1[40] -3.25 -3.24 0.91 0.86 -4.74 -1.78 1.00 2024 1897
## y1[41] 3.37 3.36 0.92 0.92 1.90 4.90 1.00 1967 1931
## y1[42] 2.41 2.41 0.90 0.90 0.93 3.89 1.00 2144 2020
## y1[43] -0.03 -0.03 0.90 0.88 -1.54 1.47 1.00 2038 2054
## y1[44] 4.32 4.33 0.92 0.94 2.84 5.82 1.00 1853 1973
## y1[45] -1.74 -1.72 0.90 0.90 -3.22 -0.31 1.00 1941 2000
## y1[46] 2.99 3.00 0.91 0.89 1.43 4.44 1.00 2021 1955
## y1[47] -0.83 -0.82 0.91 0.90 -2.31 0.68 1.00 2074 2100
## y1[48] 1.25 1.24 0.92 0.87 -0.29 2.74 1.00 1807 1931
## y1[49] 2.44 2.46 0.91 0.90 0.91 3.88 1.00 2033 1856
## y1[50] -0.57 -0.57 0.89 0.84 -2.06 0.86 1.00 1904 1745
## m1[1] 0.87 0.87 0.20 0.19 0.54 1.19 1.00 1271 963
## m1[2] 0.67 0.68 0.19 0.19 0.35 0.99 1.00 1254 994
## m1[3] 3.59 3.59 0.27 0.26 3.16 4.04 1.00 1535 1312
## m1[4] 2.26 2.26 0.26 0.26 1.85 2.69 1.00 1401 1175
## m1[5] -2.33 -2.32 0.32 0.31 -2.85 -1.81 1.00 1820 1623
## m1[6] -1.91 -1.91 0.28 0.29 -2.37 -1.45 1.00 1568 1468
## m1[7] 1.10 1.10 0.21 0.21 0.75 1.45 1.00 1378 1238
## m1[8] -1.68 -1.68 0.22 0.21 -2.03 -1.31 1.00 1445 1087
## m1[9] -0.29 -0.28 0.24 0.23 -0.69 0.11 1.00 1524 1375
## m1[10] 1.91 1.91 0.22 0.22 1.55 2.27 1.00 1429 1212
## m1[11] 1.96 1.95 0.25 0.25 1.56 2.37 1.00 1361 1022
## m1[12] 0.53 0.53 0.22 0.22 0.15 0.90 1.00 1436 1326
## m1[13] -2.11 -2.11 0.29 0.29 -2.57 -1.64 1.00 1603 1518
## m1[14] 2.74 2.74 0.25 0.24 2.33 3.15 1.00 1595 1374
## m1[15] -3.52 -3.52 0.30 0.29 -3.99 -3.03 1.00 1806 1286
## m1[16] 0.08 0.08 0.19 0.19 -0.23 0.39 1.00 1238 991
## m1[17] 2.33 2.34 0.23 0.23 1.95 2.72 1.00 1511 1333
## m1[18] -1.00 -0.99 0.27 0.26 -1.42 -0.56 1.00 1620 1431
## m1[19] 2.84 2.85 0.25 0.25 2.43 3.26 1.00 1616 1375
## m1[20] 4.04 4.04 0.26 0.26 3.61 4.46 1.00 1650 1671
## m1[21] 1.76 1.76 0.22 0.21 1.41 2.11 1.00 1400 1168
## m1[22] -2.10 -2.11 0.29 0.29 -2.56 -1.63 1.00 1602 1518
## m1[23] 2.03 2.03 0.21 0.21 1.68 2.38 1.00 1364 1500
## m1[24] 0.03 0.04 0.19 0.19 -0.28 0.35 1.00 1239 991
## m1[25] 0.95 0.95 0.20 0.20 0.62 1.27 1.00 1279 973
## m1[26] -0.54 -0.55 0.25 0.25 -0.94 -0.14 1.00 1365 1096
## m1[27] 1.67 1.67 0.21 0.21 1.33 2.02 1.00 1384 1087
## m1[28] 1.47 1.47 0.25 0.24 1.07 1.86 1.00 1431 1126
## m1[29] -1.60 -1.60 0.33 0.30 -2.17 -1.07 1.00 1776 1602
## m1[30] 3.02 3.01 0.23 0.23 2.64 3.39 1.00 1468 1578
## m1[31] 2.53 2.53 0.25 0.24 2.13 2.93 1.00 1433 1242
## m1[32] 3.83 3.83 0.31 0.31 3.32 4.35 1.00 1626 1517
## m1[33] 4.34 4.33 0.29 0.28 3.87 4.83 1.00 1634 1313
## m1[34] 2.99 2.99 0.23 0.23 2.62 3.36 1.00 1464 1584
## m1[35] 1.59 1.59 0.21 0.21 1.24 1.93 1.00 1359 1413
## m1[36] 0.32 0.33 0.23 0.22 -0.06 0.69 1.00 1454 1199
## m1[37] -0.08 -0.08 0.24 0.23 -0.48 0.30 1.00 1500 1380
## m1[38] -1.82 -1.82 0.23 0.22 -2.18 -1.44 1.00 1475 1167
## m1[39] 3.62 3.61 0.27 0.26 3.18 4.06 1.00 1537 1312
## m1[40] -3.25 -3.25 0.29 0.28 -3.71 -2.78 1.00 1755 1251
## m1[41] 3.35 3.35 0.26 0.25 2.93 3.78 1.00 1504 1260
## m1[42] 2.45 2.45 0.25 0.24 2.05 2.85 1.00 1429 1252
## m1[43] 0.00 0.00 0.24 0.24 -0.38 0.39 1.00 1315 1003
## m1[44] 4.31 4.30 0.27 0.27 3.87 4.75 1.00 1693 1697
## m1[45] -1.71 -1.71 0.22 0.21 -2.06 -1.34 1.00 1450 1087
## m1[46] 2.99 2.99 0.23 0.23 2.61 3.36 1.00 1464 1584
## m1[47] -0.82 -0.82 0.25 0.25 -1.23 -0.41 1.00 1396 1283
## m1[48] 1.23 1.23 0.24 0.24 0.84 1.62 1.00 1305 937
## m1[49] 2.43 2.44 0.24 0.24 2.04 2.82 1.00 1533 1351
## m1[50] -0.60 -0.60 0.20 0.19 -0.92 -0.28 1.00 1278 955
fn(f2)## Initial log joint probability = -8920.61
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -12.484 9.52396e-06 0.000606877 0.2695 0.8802 51
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.91 -15.57 1.88 1.66 -19.37 -13.56 1.01 762 803
## b[1] -0.19 -0.18 0.20 0.20 -0.51 0.12 1.00 934 1094
## b[2] 1.20 1.20 0.07 0.07 1.07 1.32 1.00 2107 1772
## b[3] 0.85 0.84 0.33 0.33 0.34 1.38 1.00 691 718
## b[4] 1.19 1.19 0.30 0.30 0.68 1.69 1.00 714 904
## b[5] -0.93 -0.92 0.51 0.48 -1.75 -0.09 1.00 482 550
## s 0.84 0.84 0.09 0.09 0.71 1.01 1.00 2042 1703
## y1[1] 0.73 0.73 0.87 0.85 -0.76 2.15 1.00 1842 1705
## y1[2] 0.54 0.52 0.87 0.85 -0.88 2.00 1.00 1927 1900
## y1[3] 3.28 3.30 0.90 0.89 1.79 4.74 1.00 1935 2013
## y1[4] 2.57 2.54 0.89 0.92 1.08 4.04 1.00 1522 1773
## y1[5] -2.17 -2.16 0.90 0.89 -3.65 -0.73 1.00 2206 1906
## y1[6] -1.67 -1.68 0.89 0.86 -3.14 -0.18 1.00 2072 1883
## y1[7] 1.28 1.31 0.90 0.89 -0.22 2.76 1.00 1785 1925
## y1[8] -1.88 -1.87 0.88 0.86 -3.35 -0.46 1.00 2028 1924
## y1[9] -0.10 -0.14 0.89 0.89 -1.53 1.35 1.00 1839 1816
## y1[10] 1.82 1.82 0.87 0.88 0.39 3.23 1.00 1722 2016
## y1[11] 2.23 2.22 0.90 0.89 0.70 3.70 1.00 1742 1851
## y1[12] 0.68 0.65 0.89 0.87 -0.75 2.16 1.00 2095 1973
## y1[13] -1.90 -1.89 0.91 0.91 -3.38 -0.44 1.00 1728 1819
## y1[14] 2.63 2.62 0.88 0.89 1.15 4.07 1.00 2116 2015
## y1[15] -3.72 -3.72 0.92 0.88 -5.19 -2.21 1.00 1775 1764
## y1[16] -0.06 -0.08 0.85 0.87 -1.43 1.37 1.00 1917 1846
## y1[17] 2.22 2.22 0.86 0.83 0.79 3.61 1.00 1915 1885
## y1[18] -0.87 -0.86 0.88 0.85 -2.36 0.58 1.00 1888 1891
## y1[19] 2.75 2.72 0.90 0.84 1.34 4.28 1.00 1752 1857
## y1[20] 4.26 4.25 0.90 0.91 2.77 5.72 1.00 2008 1932
## y1[21] 1.65 1.64 0.89 0.91 0.21 3.10 1.00 2061 1859
## y1[22] -1.90 -1.88 0.90 0.87 -3.37 -0.42 1.00 1629 1996
## y1[23] 2.21 2.20 0.88 0.88 0.78 3.67 1.00 2111 1920
## y1[24] -0.09 -0.07 0.88 0.86 -1.50 1.38 1.00 1744 1895
## y1[25] 0.81 0.81 0.87 0.84 -0.60 2.30 1.00 1836 1972
## y1[26] -0.30 -0.32 0.89 0.87 -1.73 1.15 1.00 1946 1968
## y1[27] 1.54 1.56 0.88 0.86 0.11 2.96 1.00 1774 1940
## y1[28] 1.12 1.13 0.91 0.88 -0.40 2.65 1.00 1810 1884
## y1[29] -2.02 -2.01 0.93 0.91 -3.56 -0.47 1.00 1721 1915
## y1[30] 3.23 3.23 0.90 0.87 1.72 4.69 1.00 1666 1566
## y1[31] 2.21 2.22 0.91 0.90 0.70 3.68 1.00 2120 1870
## y1[32] 4.16 4.16 0.91 0.92 2.70 5.61 1.00 1550 1858
## y1[33] 4.07 4.09 0.89 0.87 2.61 5.57 1.00 1905 1930
## y1[34] 3.20 3.17 0.87 0.87 1.78 4.66 1.00 1905 1869
## y1[35] 1.76 1.75 0.82 0.82 0.45 3.13 1.00 1990 1730
## y1[36] 0.50 0.48 0.90 0.93 -0.93 2.00 1.00 1776 1766
## y1[37] 0.09 0.08 0.85 0.83 -1.32 1.48 1.00 1743 1893
## y1[38] -2.00 -1.99 0.87 0.86 -3.44 -0.54 1.00 2016 2037
## y1[39] 3.30 3.31 0.89 0.89 1.80 4.81 1.00 1936 2057
## y1[40] -3.45 -3.44 0.90 0.86 -4.92 -1.97 1.00 1728 2011
## y1[41] 3.03 3.06 0.94 0.94 1.44 4.56 1.00 1773 1563
## y1[42] 2.13 2.14 0.91 0.89 0.68 3.58 1.00 1705 1784
## y1[43] 0.24 0.26 0.90 0.89 -1.25 1.69 1.00 1594 1772
## y1[44] 4.56 4.54 0.90 0.91 3.07 6.03 1.00 1780 1930
## y1[45] -1.86 -1.86 0.88 0.87 -3.30 -0.38 1.00 1967 1873
## y1[46] 3.20 3.22 0.87 0.88 1.75 4.59 1.00 1912 1899
## y1[47] -0.58 -0.56 0.88 0.87 -2.07 0.85 1.00 1881 1957
## y1[48] 1.49 1.49 0.87 0.85 0.07 2.88 1.00 1892 1898
## y1[49] 2.29 2.30 0.86 0.84 0.83 3.64 1.00 2028 1984
## y1[50] -0.77 -0.77 0.86 0.83 -2.18 0.57 1.00 1775 1808
## m1[1] 0.74 0.74 0.20 0.20 0.41 1.06 1.00 972 1173
## m1[2] 0.54 0.55 0.20 0.19 0.21 0.86 1.00 957 1165
## m1[3] 3.28 3.29 0.31 0.30 2.79 3.77 1.00 1269 1422
## m1[4] 2.55 2.54 0.30 0.29 2.08 3.06 1.00 1006 1094
## m1[5] -2.20 -2.20 0.32 0.33 -2.72 -1.68 1.00 1730 1604
## m1[6] -1.69 -1.69 0.30 0.29 -2.17 -1.19 1.00 1277 1509
## m1[7] 1.29 1.29 0.23 0.23 0.90 1.66 1.00 1192 1322
## m1[8] -1.85 -1.84 0.23 0.23 -2.24 -1.47 1.00 1068 1431
## m1[9] -0.13 -0.12 0.25 0.25 -0.53 0.29 1.00 1368 1539
## m1[10] 1.80 1.80 0.22 0.21 1.43 2.16 1.00 1109 1232
## m1[11] 2.24 2.24 0.29 0.28 1.79 2.74 1.00 992 1115
## m1[12] 0.70 0.70 0.24 0.23 0.31 1.10 1.00 1245 1371
## m1[13] -1.90 -1.89 0.30 0.30 -2.39 -1.38 1.00 1309 1488
## m1[14] 2.64 2.64 0.25 0.25 2.22 3.06 1.00 1259 1407
## m1[15] -3.72 -3.71 0.31 0.30 -4.23 -3.20 1.00 1334 1767
## m1[16] -0.06 -0.06 0.20 0.19 -0.38 0.25 1.00 931 1132
## m1[17] 2.23 2.23 0.23 0.23 1.83 2.62 1.00 1179 1248
## m1[18] -0.85 -0.84 0.27 0.28 -1.29 -0.40 1.00 1492 1623
## m1[19] 2.75 2.75 0.25 0.25 2.32 3.17 1.00 1280 1521
## m1[20] 4.28 4.27 0.29 0.28 3.79 4.76 1.00 1292 1516
## m1[21] 1.65 1.65 0.22 0.21 1.28 2.00 1.00 1085 1208
## m1[22] -1.89 -1.88 0.30 0.30 -2.38 -1.38 1.00 1308 1488
## m1[23] 2.23 2.23 0.24 0.24 1.84 2.62 1.00 1166 1312
## m1[24] -0.11 -0.10 0.20 0.19 -0.42 0.20 1.00 933 1106
## m1[25] 0.82 0.83 0.20 0.20 0.50 1.15 1.00 979 1108
## m1[26] -0.30 -0.30 0.27 0.27 -0.74 0.15 1.00 1070 1167
## m1[27] 1.56 1.56 0.21 0.20 1.19 1.91 1.00 1072 1206
## m1[28] 1.12 1.12 0.31 0.29 0.61 1.63 1.00 1163 1198
## m1[29] -2.00 -1.99 0.40 0.38 -2.66 -1.34 1.00 1277 1255
## m1[30] 3.23 3.23 0.26 0.25 2.82 3.66 1.00 1204 1358
## m1[31] 2.20 2.21 0.30 0.29 1.71 2.70 1.00 1194 1228
## m1[32] 4.14 4.14 0.35 0.35 3.58 4.73 1.00 1112 1287
## m1[33] 4.04 4.06 0.32 0.32 3.52 4.55 1.00 1356 1463
## m1[34] 3.21 3.21 0.26 0.25 2.79 3.63 1.00 1202 1360
## m1[35] 1.78 1.78 0.23 0.23 1.40 2.16 1.00 1171 1270
## m1[36] 0.50 0.50 0.24 0.23 0.10 0.89 1.00 1270 1401
## m1[37] 0.08 0.08 0.25 0.25 -0.32 0.48 1.00 1334 1474
## m1[38] -1.99 -1.99 0.24 0.23 -2.39 -1.61 1.00 1088 1459
## m1[39] 3.31 3.31 0.31 0.30 2.81 3.80 1.00 1271 1444
## m1[40] -3.45 -3.44 0.30 0.29 -3.94 -2.95 1.00 1302 1705
## m1[41] 3.04 3.04 0.30 0.30 2.54 3.53 1.00 1248 1379
## m1[42] 2.12 2.12 0.30 0.29 1.62 2.62 1.00 1191 1205
## m1[43] 0.25 0.25 0.27 0.27 -0.18 0.69 1.00 1019 1072
## m1[44] 4.54 4.54 0.30 0.28 4.05 5.04 1.00 1310 1539
## m1[45] -1.88 -1.88 0.23 0.23 -2.27 -1.50 1.00 1072 1387
## m1[46] 3.21 3.20 0.26 0.25 2.79 3.63 1.00 1202 1360
## m1[47] -0.59 -0.59 0.27 0.27 -1.02 -0.13 1.00 1105 1226
## m1[48] 1.50 1.50 0.28 0.27 1.06 1.96 1.00 978 1070
## m1[49] 2.33 2.33 0.24 0.23 1.92 2.73 1.00 1198 1303
## m1[50] -0.75 -0.75 0.20 0.20 -1.09 -0.43 1.00 948 1185
tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,xa,y,col=ca),
qplot(data=tb,xa,y,col=cb),
qplot(data=tb,xb,y,col=ca),
qplot(data=tb,xb,y,col=cb))
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)
fn(f0)## Initial log joint probability = -127.979
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 28 -49.6612 0.000110023 0.00152751 1 1 33
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -52.38 -51.97 1.89 1.72 -55.90 -50.03 1.00 799 1344
## b[1] -0.05 -0.05 0.42 0.41 -0.74 0.65 1.00 761 832
## b[2] 1.16 1.16 0.24 0.22 0.76 1.54 1.00 1996 1109
## b[3] 1.18 1.18 0.22 0.22 0.82 1.54 1.00 2208 1446
## b[4] 2.15 2.14 0.51 0.49 1.32 2.97 1.01 960 985
## b[5] 1.79 1.80 0.51 0.49 0.93 2.63 1.00 745 982
## s 1.78 1.76 0.20 0.18 1.49 2.14 1.00 1916 1374
## y1[1] 5.17 5.12 1.85 1.82 2.16 8.29 1.00 1778 1893
## y1[2] 3.92 3.91 1.86 1.83 0.95 7.00 1.00 1844 1802
## y1[3] 4.77 4.77 1.93 1.96 1.59 8.02 1.00 1824 1860
## y1[4] 3.76 3.74 1.91 1.93 0.59 6.85 1.00 1929 1974
## y1[5] 4.58 4.56 1.87 1.79 1.44 7.59 1.00 1988 1977
## y1[6] 0.09 0.07 1.85 1.84 -2.99 3.13 1.00 1947 1941
## y1[7] 5.98 5.99 1.88 1.87 2.89 9.14 1.00 1895 2012
## y1[8] 2.09 2.12 1.84 1.76 -0.95 5.14 1.00 1828 1980
## y1[9] -1.08 -1.06 1.85 1.85 -4.16 1.93 1.00 1830 2012
## y1[10] 2.30 2.34 1.84 1.79 -0.82 5.35 1.00 2023 1981
## y1[11] 4.77 4.78 1.93 1.91 1.51 7.91 1.00 1938 1936
## y1[12] 2.09 2.06 1.90 1.90 -1.00 5.24 1.00 2075 1848
## y1[13] 4.10 4.11 1.90 1.88 0.96 7.10 1.00 1776 1893
## y1[14] 4.33 4.31 1.94 1.90 1.24 7.52 1.00 1793 1742
## y1[15] 1.15 1.11 1.84 1.79 -1.77 4.21 1.00 1928 1939
## y1[16] 0.95 0.95 1.89 1.87 -2.17 4.04 1.00 2018 2054
## y1[17] 1.68 1.72 1.85 1.79 -1.46 4.68 1.00 1739 1718
## y1[18] 1.16 1.13 1.89 1.84 -1.88 4.31 1.00 1862 1842
## y1[19] -0.68 -0.68 1.90 1.85 -3.76 2.40 1.00 1957 1895
## y1[20] 4.15 4.15 1.81 1.79 1.22 7.07 1.00 1875 1746
## y1[21] 0.54 0.52 1.87 1.83 -2.52 3.57 1.00 2049 1872
## y1[22] 3.01 3.01 1.86 1.82 0.10 6.12 1.00 1709 1704
## y1[23] -0.48 -0.42 1.81 1.76 -3.62 2.41 1.00 1816 1917
## y1[24] 3.00 2.97 1.83 1.85 0.03 5.82 1.00 2058 1811
## y1[25] 4.21 4.22 1.85 1.86 1.24 7.23 1.00 2073 1918
## y1[26] 2.00 1.99 1.87 1.88 -1.00 5.14 1.00 1918 1864
## y1[27] 1.06 1.05 1.90 1.86 -2.05 4.13 1.00 1860 1666
## y1[28] 0.17 0.18 1.86 1.93 -2.90 3.17 1.00 1739 1841
## y1[29] -1.10 -1.06 1.87 1.83 -4.12 1.91 1.00 1838 1946
## y1[30] -1.92 -1.92 1.92 1.87 -5.01 1.21 1.00 1751 1790
## y1[31] 2.30 2.35 1.87 1.86 -0.85 5.43 1.00 1956 1864
## y1[32] -0.94 -0.93 1.86 1.80 -4.04 2.10 1.00 2010 1895
## y1[33] -1.87 -1.87 1.89 1.92 -4.94 1.19 1.00 1954 1744
## y1[34] -0.61 -0.60 1.86 1.81 -3.65 2.48 1.00 1788 1665
## y1[35] 2.26 2.24 1.89 1.99 -0.80 5.30 1.00 2035 1917
## y1[36] -0.73 -0.73 1.91 1.88 -3.87 2.42 1.00 1895 1974
## y1[37] -2.15 -2.13 1.95 1.88 -5.41 1.05 1.00 1890 1857
## y1[38] 3.60 3.58 1.83 1.78 0.67 6.67 1.00 2116 2044
## y1[39] 1.29 1.34 1.84 1.84 -1.81 4.34 1.00 1785 1893
## y1[40] 0.09 0.05 1.85 1.89 -2.84 3.12 1.00 1940 1971
## y1[41] 0.35 0.30 1.90 1.88 -2.76 3.48 1.00 2054 1973
## y1[42] 3.65 3.60 1.83 1.75 0.60 6.70 1.00 1956 2015
## y1[43] 3.64 3.65 1.88 1.83 0.46 6.61 1.00 1955 1805
## y1[44] 0.78 0.80 1.81 1.78 -2.21 3.77 1.00 1900 2014
## y1[45] 1.05 1.08 1.89 1.84 -2.10 4.12 1.00 1898 1856
## y1[46] 2.02 2.05 1.84 1.79 -0.99 5.10 1.00 1855 1892
## y1[47] -1.24 -1.20 1.91 1.95 -4.28 1.79 1.00 1706 1823
## y1[48] 1.30 1.30 1.85 1.86 -1.76 4.42 1.00 1813 1967
## y1[49] 3.95 3.96 1.89 1.87 0.78 7.04 1.00 1846 1796
## y1[50] 1.25 1.28 1.86 1.84 -1.71 4.19 1.00 1644 2014
## m1[1] 5.18 5.18 0.55 0.53 4.26 6.07 1.01 1460 1528
## m1[2] 3.86 3.85 0.55 0.52 2.93 4.74 1.00 1598 1548
## m1[3] 4.79 4.77 0.61 0.60 3.78 5.79 1.00 1473 1471
## m1[4] 3.66 3.65 0.61 0.59 2.64 4.67 1.00 1537 1406
## m1[5] 4.51 4.51 0.49 0.48 3.69 5.30 1.00 1407 1726
## m1[6] 0.14 0.14 0.53 0.53 -0.74 1.02 1.00 1822 1524
## m1[7] 5.93 5.93 0.64 0.61 4.85 6.96 1.00 1541 1362
## m1[8] 2.14 2.14 0.52 0.51 1.28 3.01 1.00 1034 1166
## m1[9] -1.07 -1.07 0.59 0.59 -2.08 -0.11 1.00 1907 1719
## m1[10] 2.27 2.26 0.49 0.49 1.51 3.09 1.00 984 1244
## m1[11] 4.80 4.80 0.62 0.60 3.79 5.80 1.00 2135 1539
## m1[12] 2.06 2.08 0.72 0.69 0.90 3.28 1.00 1934 1542
## m1[13] 4.14 4.15 0.67 0.64 3.00 5.21 1.00 1528 1430
## m1[14] 4.30 4.33 0.69 0.66 3.19 5.41 1.00 1707 1463
## m1[15] 1.22 1.21 0.50 0.48 0.39 2.06 1.00 954 1178
## m1[16] 0.92 0.91 0.55 0.54 0.02 1.84 1.00 1953 1362
## m1[17] 1.67 1.68 0.54 0.53 0.77 2.57 1.00 1369 1417
## m1[18] 1.14 1.15 0.65 0.62 0.07 2.21 1.00 1341 1340
## m1[19] -0.67 -0.67 0.59 0.58 -1.64 0.27 1.00 1610 1551
## m1[20] 4.14 4.15 0.54 0.52 3.23 4.99 1.00 1409 1313
## m1[21] 0.55 0.55 0.61 0.61 -0.45 1.56 1.00 1402 1389
## m1[22] 3.00 2.98 0.57 0.57 2.09 3.93 1.00 1158 1321
## m1[23] -0.48 -0.48 0.54 0.52 -1.39 0.39 1.00 1476 1355
## m1[24] 3.01 3.00 0.48 0.46 2.23 3.80 1.00 1275 1412
## m1[25] 4.26 4.25 0.48 0.47 3.47 5.03 1.00 1345 1311
## m1[26] 1.99 1.97 0.51 0.49 1.17 2.85 1.00 1654 1290
## m1[27] 1.05 1.06 0.56 0.55 0.13 1.98 1.00 1818 1634
## m1[28] 0.18 0.18 0.51 0.51 -0.68 1.02 1.00 1526 1361
## m1[29] -1.11 -1.12 0.62 0.62 -2.12 -0.08 1.00 1199 1362
## m1[30] -1.93 -1.92 0.69 0.65 -3.08 -0.84 1.00 1858 1629
## m1[31] 2.26 2.26 0.57 0.56 1.32 3.21 1.00 1406 1385
## m1[32] -0.98 -0.98 0.59 0.58 -1.99 -0.04 1.00 1890 1688
## m1[33] -1.84 -1.84 0.55 0.53 -2.74 -0.94 1.00 1035 1342
## m1[34] -0.67 -0.67 0.61 0.60 -1.67 0.32 1.00 1653 1511
## m1[35] 2.27 2.26 0.49 0.49 1.50 3.08 1.00 981 1181
## m1[36] -0.68 -0.68 0.68 0.68 -1.80 0.42 1.00 1144 1396
## m1[37] -2.12 -2.11 0.69 0.68 -3.28 -1.01 1.00 1584 1304
## m1[38] 3.65 3.65 0.50 0.49 2.85 4.48 1.00 1784 1694
## m1[39] 1.35 1.34 0.48 0.47 0.57 2.14 1.00 909 1166
## m1[40] 0.10 0.10 0.50 0.49 -0.75 0.94 1.00 1589 1571
## m1[41] 0.32 0.32 0.66 0.65 -0.74 1.45 1.00 1316 1257
## m1[42] 3.64 3.63 0.46 0.44 2.89 4.39 1.00 1369 1187
## m1[43] 3.67 3.66 0.53 0.51 2.81 4.56 1.00 1398 1385
## m1[44] 0.79 0.80 0.47 0.46 0.01 1.58 1.00 1439 1431
## m1[45] 1.05 1.05 0.58 0.56 0.10 2.01 1.00 1055 1183
## m1[46] 2.06 2.06 0.45 0.44 1.32 2.79 1.00 1557 1566
## m1[47] -1.27 -1.26 0.61 0.59 -2.27 -0.28 1.00 1515 1385
## m1[48] 1.33 1.33 0.43 0.44 0.64 2.06 1.00 830 1031
## m1[49] 3.95 3.94 0.55 0.53 3.06 4.87 1.00 1626 1520
## m1[50] 1.33 1.32 0.48 0.47 0.55 2.12 1.00 908 1180
fn(f1)## Initial log joint probability = -125.882
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 35 -49.6599 0.000110695 0.0017254 0.4031 1 41
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -53.07 -52.70 2.13 1.94 -56.99 -50.29 1.00 544 945
## b[1] -0.03 -0.03 0.54 0.54 -0.92 0.88 1.00 687 880
## b[2] 1.16 1.16 0.24 0.24 0.76 1.56 1.00 2298 1450
## b[3] 1.17 1.17 0.24 0.24 0.78 1.56 1.00 1659 1536
## b[4] 2.07 2.07 0.82 0.80 0.74 3.40 1.01 581 778
## b[5] 1.77 1.78 0.76 0.75 0.51 3.01 1.01 640 847
## b[6] 0.07 0.02 1.17 1.06 -1.83 1.93 1.01 548 765
## s 1.80 1.79 0.20 0.20 1.51 2.16 1.00 2078 1346
## y1[1] 5.09 5.11 1.94 1.97 1.92 8.23 1.00 1975 1806
## y1[2] 3.84 3.87 1.93 1.85 0.61 6.87 1.00 1934 1969
## y1[3] 4.75 4.75 1.97 1.98 1.42 7.97 1.00 1734 1700
## y1[4] 3.71 3.67 1.96 1.98 0.53 6.82 1.00 1962 1916
## y1[5] 4.50 4.48 1.88 1.86 1.50 7.57 1.00 2105 1840
## y1[6] 0.07 0.06 1.90 1.85 -3.15 3.22 1.00 1804 1909
## y1[7] 5.97 6.00 1.92 1.92 2.77 9.15 1.00 2062 2015
## y1[8] 2.23 2.28 1.85 1.83 -0.83 5.23 1.00 2169 1990
## y1[9] -1.09 -1.10 1.91 1.94 -4.23 1.95 1.00 1893 1794
## y1[10] 2.32 2.31 1.88 1.87 -0.76 5.43 1.00 1961 2009
## y1[11] 4.82 4.77 2.03 2.06 1.58 8.18 1.00 2007 1973
## y1[12] 2.03 2.02 1.95 1.85 -1.29 5.22 1.00 2172 1915
## y1[13] 4.07 4.07 1.99 1.83 0.86 7.43 1.00 1998 2088
## y1[14] 4.21 4.16 1.92 1.88 1.10 7.38 1.00 2106 1894
## y1[15] 1.24 1.20 1.91 1.86 -1.88 4.35 1.00 1657 1804
## y1[16] 0.87 0.85 1.95 1.95 -2.33 4.12 1.00 1875 1885
## y1[17] 1.68 1.71 1.92 1.96 -1.58 4.71 1.00 1849 2000
## y1[18] 1.09 1.14 1.92 1.87 -2.05 4.26 1.00 1977 1855
## y1[19] -0.69 -0.65 1.91 1.95 -3.78 2.43 1.00 1972 1807
## y1[20] 4.12 4.08 1.91 1.89 0.97 7.28 1.00 2057 1916
## y1[21] 0.56 0.60 1.90 1.80 -2.63 3.63 1.00 1755 1909
## y1[22] 3.04 3.05 1.89 1.83 -0.05 6.19 1.00 1966 1919
## y1[23] -0.54 -0.54 1.89 1.88 -3.64 2.43 1.00 2144 1836
## y1[24] 2.92 2.92 1.93 1.93 -0.25 6.01 1.00 1722 1820
## y1[25] 4.25 4.29 1.89 1.83 1.08 7.35 1.00 1985 1820
## y1[26] 1.94 2.00 1.95 1.86 -1.33 5.14 1.00 1919 1838
## y1[27] 1.03 1.03 1.93 1.97 -2.14 4.15 1.00 1741 1917
## y1[28] 0.15 0.17 1.87 1.89 -2.96 3.16 1.00 1770 1788
## y1[29] -1.04 -0.99 1.90 1.81 -4.18 1.99 1.00 1714 1725
## y1[30] -1.96 -1.92 1.97 1.96 -5.12 1.26 1.00 2026 1815
## y1[31] 2.30 2.18 1.97 1.94 -0.82 5.74 1.00 2011 1922
## y1[32] -1.03 -0.92 1.96 1.94 -4.25 2.06 1.00 1910 1836
## y1[33] -1.84 -1.85 1.94 1.92 -4.96 1.49 1.00 1500 1682
## y1[34] -0.73 -0.73 1.89 1.80 -3.98 2.34 1.00 2073 1949
## y1[35] 2.25 2.20 1.87 1.85 -0.76 5.30 1.00 1809 1762
## y1[36] -0.67 -0.69 1.98 1.98 -3.89 2.63 1.00 2035 1969
## y1[37] -2.11 -2.07 1.93 1.91 -5.34 0.93 1.00 2126 1901
## y1[38] 3.56 3.55 1.92 1.86 0.38 6.67 1.00 1912 1894
## y1[39] 1.31 1.35 1.89 1.89 -1.75 4.37 1.00 1890 1782
## y1[40] 0.09 0.12 1.86 1.81 -2.97 3.11 1.00 1785 1985
## y1[41] 0.37 0.41 1.97 1.93 -2.95 3.60 1.00 1560 1861
## y1[42] 3.59 3.59 1.85 1.81 0.50 6.58 1.00 1994 2057
## y1[43] 3.62 3.62 1.93 1.92 0.42 6.72 1.00 1660 1840
## y1[44] 0.74 0.75 1.88 1.89 -2.41 3.74 1.00 2134 2009
## y1[45] 1.10 1.10 1.98 1.89 -2.22 4.28 1.00 1960 2016
## y1[46] 2.07 2.07 1.88 1.89 -0.97 5.21 1.00 1886 1875
## y1[47] -1.22 -1.20 1.91 1.91 -4.36 1.96 1.00 2050 1910
## y1[48] 1.38 1.35 1.92 1.88 -1.75 4.57 1.00 1983 1809
## y1[49] 3.93 3.94 1.90 1.87 0.82 7.09 1.00 1879 1800
## y1[50] 1.36 1.37 1.87 1.83 -1.74 4.44 1.00 1643 1854
## m1[1] 5.16 5.16 0.62 0.61 4.16 6.19 1.00 2063 1256
## m1[2] 3.84 3.85 0.64 0.61 2.80 4.86 1.00 2045 1196
## m1[3] 4.73 4.72 0.76 0.73 3.46 5.97 1.00 979 1167
## m1[4] 3.67 3.67 0.69 0.66 2.53 4.78 1.00 1549 1480
## m1[5] 4.50 4.49 0.58 0.56 3.55 5.44 1.00 1904 1660
## m1[6] 0.10 0.11 0.55 0.55 -0.79 1.01 1.00 1582 1482
## m1[7] 5.91 5.91 0.70 0.68 4.79 7.05 1.00 2200 1220
## m1[8] 2.15 2.14 0.57 0.55 1.22 3.12 1.00 1296 1542
## m1[9] -1.11 -1.12 0.60 0.60 -2.12 -0.13 1.00 1863 1664
## m1[10] 2.29 2.28 0.54 0.53 1.37 3.18 1.00 1164 1303
## m1[11] 4.78 4.78 0.78 0.79 3.52 6.01 1.00 1137 1354
## m1[12] 2.03 2.02 0.74 0.71 0.81 3.25 1.00 1715 1375
## m1[13] 4.14 4.12 0.77 0.76 2.87 5.41 1.00 1729 1857
## m1[14] 4.27 4.27 0.76 0.73 3.03 5.52 1.00 2317 1207
## m1[15] 1.23 1.22 0.56 0.55 0.27 2.19 1.00 1057 1384
## m1[16] 0.92 0.92 0.73 0.69 -0.30 2.12 1.00 1244 1415
## m1[17] 1.68 1.67 0.58 0.57 0.73 2.63 1.00 1820 1615
## m1[18] 1.07 1.06 0.74 0.72 -0.12 2.27 1.00 1490 1327
## m1[19] -0.68 -0.67 0.61 0.60 -1.77 0.31 1.00 1868 1625
## m1[20] 4.13 4.12 0.65 0.65 3.08 5.19 1.00 1580 1723
## m1[21] 0.57 0.57 0.62 0.61 -0.45 1.59 1.00 2019 1681
## m1[22] 3.02 3.02 0.61 0.58 1.98 4.04 1.00 1520 1376
## m1[23] -0.48 -0.46 0.55 0.54 -1.43 0.42 1.00 1914 1382
## m1[24] 2.96 2.97 0.60 0.57 1.97 3.91 1.00 950 1075
## m1[25] 4.25 4.24 0.58 0.59 3.32 5.16 1.00 1632 1634
## m1[26] 1.99 1.99 0.67 0.64 0.86 3.07 1.00 1272 1494
## m1[27] 1.01 1.03 0.59 0.58 0.08 1.96 1.00 1528 1146
## m1[28] 0.18 0.18 0.55 0.54 -0.79 1.06 1.00 1632 1630
## m1[29] -1.08 -1.09 0.78 0.74 -2.37 0.22 1.00 842 1067
## m1[30] -1.98 -1.98 0.70 0.70 -3.11 -0.83 1.00 2096 1943
## m1[31] 2.26 2.26 0.61 0.61 1.25 3.29 1.00 1819 1467
## m1[32] -1.03 -1.03 0.60 0.60 -2.02 -0.06 1.00 1836 1664
## m1[33] -1.81 -1.82 0.72 0.70 -2.99 -0.63 1.00 696 1009
## m1[34] -0.68 -0.66 0.64 0.62 -1.81 0.35 1.00 1855 1666
## m1[35] 2.28 2.28 0.54 0.53 1.38 3.18 1.00 1177 1397
## m1[36] -0.68 -0.68 0.75 0.71 -1.96 0.61 1.00 1114 1382
## m1[37] -2.11 -2.10 0.69 0.67 -3.28 -1.01 1.00 2044 1572
## m1[38] 3.64 3.66 0.63 0.63 2.61 4.63 1.00 1139 1320
## m1[39] 1.36 1.35 0.54 0.53 0.44 2.28 1.00 1003 1423
## m1[40] 0.05 0.06 0.54 0.52 -0.85 0.91 1.00 1464 1437
## m1[41] 0.36 0.37 0.77 0.74 -0.94 1.64 1.00 1056 1313
## m1[42] 3.63 3.62 0.58 0.58 2.68 4.55 1.00 1480 1641
## m1[43] 3.62 3.63 0.65 0.60 2.52 4.71 1.00 976 1104
## m1[44] 0.74 0.76 0.52 0.49 -0.13 1.57 1.00 1250 1397
## m1[45] 1.06 1.04 0.63 0.62 0.01 2.13 1.00 1244 1447
## m1[46] 2.05 2.06 0.54 0.54 1.15 2.90 1.00 1235 1411
## m1[47] -1.26 -1.24 0.61 0.58 -2.30 -0.29 1.00 2048 1600
## m1[48] 1.35 1.34 0.51 0.51 0.49 2.21 1.00 862 1044
## m1[49] 3.95 3.97 0.66 0.64 2.86 5.03 1.00 1279 1462
## m1[50] 1.34 1.33 0.54 0.53 0.42 2.26 1.00 999 1423
fn(f2)## Initial log joint probability = -101.671
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -12.1861 6.24919e-05 0.000665711 0.796 0.796 32
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.76 -16.46 2.09 1.97 -20.59 -14.01 1.00 831 1308
## b[1] -0.16 -0.17 0.25 0.25 -0.57 0.23 1.00 623 840
## b[2] 0.68 0.68 0.12 0.12 0.47 0.87 1.00 2306 1607
## b[3] 1.28 1.29 0.11 0.11 1.10 1.46 1.00 1312 1545
## b[4] 2.58 2.60 0.38 0.37 1.92 3.19 1.01 513 743
## b[5] 2.39 2.40 0.36 0.35 1.77 2.96 1.01 592 717
## b[6] -1.24 -1.24 0.10 0.11 -1.41 -1.07 1.00 2235 1707
## b[7] -1.39 -1.41 0.53 0.53 -2.26 -0.49 1.01 477 695
## s 0.86 0.85 0.09 0.09 0.72 1.02 1.00 1744 1386
## y1[1] 4.82 4.83 0.89 0.87 3.36 6.27 1.00 1937 1931
## y1[2] 4.84 4.80 0.92 0.90 3.37 6.34 1.00 1950 1919
## y1[3] 3.01 3.01 0.97 0.96 1.46 4.60 1.00 1725 1958
## y1[4] 3.62 3.61 0.91 0.87 2.10 5.12 1.00 1995 1766
## y1[5] 4.10 4.09 0.88 0.87 2.66 5.54 1.00 1896 2058
## y1[6] 0.29 0.29 0.90 0.85 -1.17 1.79 1.00 1699 1657
## y1[7] 5.55 5.56 0.91 0.92 4.10 7.06 1.00 2084 1946
## y1[8] 1.32 1.33 0.92 0.89 -0.26 2.81 1.00 2168 1940
## y1[9] -2.54 -2.52 0.90 0.90 -4.01 -1.11 1.00 2111 1733
## y1[10] 0.49 0.50 0.91 0.89 -0.99 2.00 1.00 1735 1793
## y1[11] 3.63 3.64 0.95 0.91 2.07 5.14 1.00 1749 1806
## y1[12] 5.72 5.70 0.98 0.97 4.14 7.35 1.00 1829 1874
## y1[13] 5.50 5.50 0.95 0.94 3.96 7.03 1.00 1992 1934
## y1[14] 7.29 7.29 0.97 0.96 5.70 8.83 1.00 1906 1707
## y1[15] 1.53 1.53 0.89 0.88 0.10 2.98 1.00 1843 1973
## y1[16] -1.17 -1.16 0.93 0.94 -2.69 0.34 1.00 1867 1947
## y1[17] 2.73 2.73 0.91 0.95 1.24 4.19 1.00 1678 1875
## y1[18] 3.70 3.69 0.98 0.97 2.14 5.31 1.00 1621 1545
## y1[19] -0.37 -0.38 0.90 0.91 -1.80 1.07 1.00 2079 1930
## y1[20] 4.05 4.04 0.91 0.89 2.56 5.54 1.00 1959 1896
## y1[21] 1.93 1.91 0.92 0.92 0.47 3.49 1.00 1957 1718
## y1[22] 0.32 0.32 0.94 0.95 -1.20 1.85 1.00 1841 2017
## y1[23] -0.67 -0.69 0.90 0.88 -2.19 0.81 1.00 1914 1746
## y1[24] 2.95 2.95 0.92 0.90 1.45 4.45 1.00 1920 1931
## y1[25] 3.62 3.61 0.91 0.90 2.18 5.15 1.00 1804 1997
## y1[26] 1.11 1.10 0.94 0.90 -0.45 2.64 1.00 1889 1870
## y1[27] 2.38 2.38 0.88 0.89 0.93 3.82 1.00 2101 2039
## y1[28] 0.91 0.91 0.92 0.93 -0.66 2.41 1.00 2095 1704
## y1[29] 0.39 0.37 0.96 0.90 -1.19 2.01 1.00 1610 1675
## y1[30] -4.53 -4.54 0.95 0.93 -6.09 -2.92 1.00 1853 1791
## y1[31] 3.38 3.39 0.92 0.91 1.78 4.87 1.00 1800 1856
## y1[32] -2.28 -2.27 0.90 0.86 -3.78 -0.81 1.00 1971 1859
## y1[33] -2.12 -2.13 0.91 0.89 -3.67 -0.64 1.00 1848 1780
## y1[34] 0.07 0.10 0.92 0.91 -1.47 1.55 1.00 1621 1702
## y1[35] 0.56 0.57 0.91 0.92 -0.94 2.08 1.00 1858 1845
## y1[36] 2.51 2.51 0.98 0.97 0.88 4.04 1.00 2032 2015
## y1[37] -4.48 -4.46 0.97 0.98 -6.11 -2.91 1.00 1858 1737
## y1[38] 3.29 3.33 0.92 0.87 1.77 4.76 1.00 1830 1949
## y1[39] 1.19 1.19 0.87 0.90 -0.21 2.60 1.00 1552 1817
## y1[40] -0.22 -0.22 0.88 0.86 -1.67 1.23 1.00 1838 1844
## y1[41] 2.57 2.59 0.94 0.91 0.99 4.10 1.00 1608 1879
## y1[42] 3.11 3.07 0.90 0.90 1.67 4.59 1.00 1871 1820
## y1[43] 3.21 3.17 0.92 0.91 1.72 4.74 1.00 1896 2059
## y1[44] 0.87 0.89 0.87 0.84 -0.56 2.33 1.00 1954 1876
## y1[45] 2.63 2.64 0.89 0.88 1.20 4.13 1.00 1803 1972
## y1[46] 2.90 2.92 0.92 0.90 1.37 4.37 1.00 2061 1875
## y1[47] -2.40 -2.39 0.90 0.90 -3.86 -0.91 1.00 2064 2037
## y1[48] 0.50 0.53 0.90 0.88 -1.00 1.95 1.00 1669 1821
## y1[49] 2.90 2.90 0.91 0.92 1.45 4.40 1.00 1753 1858
## y1[50] 1.19 1.20 0.91 0.92 -0.31 2.65 1.00 1705 1808
## m1[1] 4.81 4.81 0.28 0.29 4.36 5.28 1.00 2095 1727
## m1[2] 4.82 4.82 0.30 0.30 4.33 5.31 1.00 2252 1571
## m1[3] 3.03 3.04 0.39 0.39 2.39 3.67 1.00 1258 1371
## m1[4] 3.59 3.59 0.32 0.32 3.08 4.13 1.00 1254 1281
## m1[5] 4.10 4.10 0.26 0.27 3.66 4.53 1.00 2022 1553
## m1[6] 0.29 0.29 0.26 0.25 -0.14 0.72 1.00 1719 1833
## m1[7] 5.55 5.55 0.31 0.32 5.02 6.04 1.00 2180 1502
## m1[8] 1.33 1.32 0.27 0.27 0.89 1.77 1.00 1106 1505
## m1[9] -2.53 -2.52 0.31 0.31 -3.04 -2.00 1.00 2145 1787
## m1[10] 0.52 0.51 0.29 0.29 0.06 1.00 1.00 1125 1444
## m1[11] 3.64 3.64 0.37 0.38 3.02 4.25 1.00 1242 1370
## m1[12] 5.73 5.73 0.48 0.47 4.98 6.49 1.00 1694 1540
## m1[13] 5.52 5.52 0.38 0.38 4.90 6.15 1.00 1933 1589
## m1[14] 7.26 7.26 0.42 0.42 6.57 7.93 1.00 2524 1499
## m1[15] 1.56 1.56 0.26 0.27 1.14 1.98 1.00 953 1209
## m1[16] -1.13 -1.13 0.38 0.39 -1.78 -0.50 1.00 1212 1632
## m1[17] 2.74 2.74 0.29 0.29 2.27 3.22 1.00 1227 1368
## m1[18] 3.69 3.69 0.40 0.39 3.03 4.34 1.00 1139 1321
## m1[19] -0.35 -0.35 0.30 0.30 -0.84 0.15 1.00 1810 1562
## m1[20] 4.04 4.04 0.29 0.30 3.54 4.52 1.00 1700 1557
## m1[21] 1.93 1.94 0.32 0.31 1.41 2.45 1.00 1383 1494
## m1[22] 0.33 0.32 0.36 0.37 -0.24 0.92 1.00 1396 1385
## m1[23] -0.66 -0.66 0.28 0.27 -1.12 -0.20 1.00 1646 1649
## m1[24] 2.92 2.93 0.28 0.27 2.46 3.39 1.00 949 1191
## m1[25] 3.63 3.63 0.27 0.27 3.18 4.06 1.00 1749 1506
## m1[26] 1.10 1.10 0.32 0.31 0.57 1.61 1.00 1292 1530
## m1[27] 2.38 2.38 0.30 0.28 1.89 2.87 1.00 1514 1566
## m1[28] 0.90 0.91 0.28 0.27 0.46 1.34 1.00 1426 1527
## m1[29] 0.38 0.38 0.38 0.35 -0.24 1.01 1.00 879 1100
## m1[30] -4.55 -4.54 0.41 0.41 -5.22 -3.85 1.00 2252 1911
## m1[31] 3.36 3.36 0.31 0.31 2.86 3.87 1.00 1224 1368
## m1[32] -2.33 -2.32 0.31 0.30 -2.83 -1.81 1.00 2117 1823
## m1[33] -2.13 -2.13 0.33 0.33 -2.65 -1.58 1.00 620 995
## m1[34] 0.06 0.07 0.32 0.32 -0.47 0.59 1.00 1818 1650
## m1[35] 0.58 0.58 0.28 0.28 0.13 1.06 1.00 1120 1401
## m1[36] 2.51 2.50 0.43 0.43 1.80 3.21 1.00 1816 1561
## m1[37] -4.46 -4.46 0.41 0.39 -5.13 -3.78 1.00 2024 1493
## m1[38] 3.26 3.25 0.30 0.30 2.77 3.74 1.00 1085 1345
## m1[39] 1.18 1.17 0.25 0.26 0.78 1.59 1.00 866 1401
## m1[40] -0.23 -0.22 0.26 0.25 -0.66 0.18 1.00 1614 1591
## m1[41] 2.58 2.57 0.41 0.40 1.89 3.26 1.00 1268 1462
## m1[42] 3.12 3.13 0.27 0.26 2.68 3.55 1.00 1608 1611
## m1[43] 3.22 3.23 0.31 0.30 2.71 3.72 1.00 1052 1304
## m1[44] 0.87 0.88 0.25 0.24 0.47 1.27 1.00 1294 1520
## m1[45] 2.62 2.62 0.31 0.33 2.11 3.15 1.00 1502 1379
## m1[46] 2.89 2.89 0.26 0.26 2.46 3.31 1.00 1009 1121
## m1[47] -2.39 -2.39 0.32 0.31 -2.93 -1.85 1.00 1902 1572
## m1[48] 0.52 0.51 0.24 0.25 0.13 0.92 1.00 779 1186
## m1[49] 2.95 2.95 0.32 0.32 2.44 3.48 1.00 1288 1441
## m1[50] 1.19 1.18 0.25 0.26 0.79 1.59 1.00 863 1401
generalized linear regression
log normal regression
objective variable [0,Infinity)
# of samples is N,
log mi=sum(bj*xji),j=[0,K],i=[1,N]
log yi~N(mi,s) ex5-2.stan
log normal regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~lognormal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,log(x1+x2),1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex5-2.stan')
mle=mdl$optimize(data=data) # it sometimes occur error and stop process## Initial log joint probability = -3273.07
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 36 -19.7607 8.31624e-05 0.00157544 1 1 57
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -19.76
## b[1] -1.22
## b[2] 0.23
## b[3] 0.15
## s 0.65
## y1[1] 0.92
## y1[2] 2.79
## y1[3] 0.55
## y1[4] 2.40
## y1[5] 2.62
## y1[6] 0.90
## y1[7] 1.96
## y1[8] 0.11
## y1[9] 3.75
## y1[10] 0.95
## y1[11] 1.92
## y1[12] 2.02
## y1[13] 5.38
## y1[14] 1.31
## y1[15] 0.80
## y1[16] 0.92
## y1[17] 0.65
## y1[18] 2.13
## y1[19] 1.81
## y1[20] 2.48
## m1[1] -0.74
## m1[2] 1.00
## m1[3] -0.12
## m1[4] 1.20
## m1[5] 0.68
## m1[6] 0.21
## m1[7] 0.56
## m1[8] -0.76
## m1[9] 1.94
## m1[10] 0.58
## m1[11] 0.33
## m1[12] 0.40
## m1[13] 1.04
## m1[14] 1.23
## m1[15] 1.28
## m1[16] 0.62
## m1[17] -0.07
## m1[18] 1.22
## m1[19] 0.49
## m1[20] 0.20mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.34 -21.98 1.51 1.36 -25.37 -20.57 1.00 555 846
## b[1] -1.24 -1.24 0.49 0.46 -2.04 -0.43 1.00 312 371
## b[2] 0.23 0.23 0.07 0.07 0.11 0.35 1.00 469 630
## b[3] 0.16 0.16 0.07 0.06 0.04 0.28 1.00 971 841
## s 0.76 0.74 0.14 0.13 0.57 1.01 1.00 861 1077
## y1[1] 0.72 0.47 0.81 0.36 0.11 2.11 1.00 1074 1196
## y1[2] 3.96 2.72 4.49 2.06 0.73 11.14 1.00 1936 1821
## y1[3] 1.25 0.88 1.29 0.63 0.23 3.32 1.00 1362 1697
## y1[4] 4.68 3.35 5.29 2.41 0.88 12.81 1.00 1867 1778
## y1[5] 2.68 1.99 2.52 1.42 0.53 7.25 1.00 2324 2098
## y1[6] 1.77 1.22 2.10 0.89 0.31 4.81 1.00 1727 1918
## y1[7] 2.42 1.76 2.42 1.24 0.50 6.25 1.00 1924 1834
## y1[8] 0.65 0.45 0.72 0.33 0.11 1.69 1.00 1032 1660
## y1[9] 10.86 7.48 12.96 5.68 1.68 31.40 1.00 1677 1772
## y1[10] 2.41 1.70 2.64 1.22 0.49 6.27 1.00 1912 1599
## y1[11] 1.98 1.42 2.05 1.10 0.36 5.43 1.00 1719 1799
## y1[12] 2.09 1.54 2.19 1.08 0.42 5.48 1.00 1664 1843
## y1[13] 4.09 2.98 4.28 2.09 0.76 11.04 1.00 1960 1931
## y1[14] 5.22 3.61 5.69 2.73 0.89 14.19 1.00 1887 2057
## y1[15] 5.10 3.63 5.38 2.63 0.96 14.00 1.00 2221 1890
## y1[16] 2.67 1.95 2.87 1.45 0.50 7.00 1.00 2047 1968
## y1[17] 1.37 0.92 1.60 0.73 0.24 3.80 1.00 1519 1682
## y1[18] 4.76 3.44 4.60 2.60 0.86 12.84 1.00 1741 1785
## y1[19] 2.47 1.64 3.19 1.25 0.39 7.07 1.00 2039 1859
## y1[20] 1.72 1.24 1.99 0.90 0.32 4.67 1.00 1933 1882
## m1[1] -0.75 -0.76 0.38 0.35 -1.40 -0.12 1.00 321 382
## m1[2] 1.01 1.01 0.27 0.26 0.59 1.45 1.00 1942 1462
## m1[3] -0.12 -0.12 0.25 0.23 -0.54 0.29 1.00 374 492
## m1[4] 1.21 1.22 0.28 0.27 0.75 1.67 1.00 1342 1382
## m1[5] 0.68 0.69 0.26 0.26 0.24 1.11 1.00 2059 1694
## m1[6] 0.21 0.22 0.29 0.28 -0.26 0.69 1.00 1088 959
## m1[7] 0.57 0.57 0.17 0.17 0.29 0.86 1.00 1396 1215
## m1[8] -0.77 -0.77 0.39 0.36 -1.42 -0.13 1.00 321 391
## m1[9] 1.96 1.96 0.42 0.41 1.28 2.64 1.00 730 958
## m1[10] 0.59 0.59 0.17 0.17 0.31 0.87 1.00 1534 1208
## m1[11] 0.33 0.33 0.28 0.27 -0.13 0.78 1.00 1377 1193
## m1[12] 0.41 0.41 0.18 0.17 0.12 0.70 1.00 803 822
## m1[13] 1.05 1.05 0.21 0.22 0.71 1.40 1.00 1938 1461
## m1[14] 1.25 1.24 0.30 0.30 0.77 1.75 1.00 1757 1493
## m1[15] 1.29 1.29 0.25 0.25 0.88 1.71 1.00 1264 1539
## m1[16] 0.63 0.63 0.32 0.31 0.09 1.13 1.00 1983 1333
## m1[17] -0.06 -0.07 0.37 0.35 -0.68 0.52 1.00 687 791
## m1[18] 1.23 1.23 0.30 0.29 0.72 1.73 1.00 1367 1475
## m1[19] 0.50 0.50 0.35 0.34 -0.08 1.08 1.00 1521 1197
## m1[20] 0.20 0.20 0.23 0.22 -0.18 0.58 1.00 689 805
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
poisson regression
objective variable [0,Infinity), integer. variance of error is near to mean
(normal linear regression, correlation is near to 1,-1, variance of error is constant)
# of samples is N,
log li=sum(bj*xji),j=[0,K],i=[1,N]
yi~Po(li)
or
li=sum(bj*xji),j=[0,k]
yi~Po(exp li)
when xj -> xj +1, y -> y* exp bj ex6-1.stan
poisson regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] l=X*b;
y~poisson_log(l);
}
generated quantities{
array[N] int y1;
vector[N] l1=X*b;
for(i in 1:N){
y1[i]=poisson_log_rng(l1[i]);
}
}n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='poisson')##
## Call: glm(formula = f, family = "poisson", data = tb)
##
## Coefficients:
## (Intercept) x cb
## -0.223 1.635 0.547
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 59.9
## Residual Deviance: 22 AIC: 81.5X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -242.607
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -4.98764 5.28446e-05 0.000456263 0.8447 0.8447 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -4.99
## b[1] -0.22
## b[2] 1.64
## b[3] 0.55
## y1[1] 1.00
## y1[2] 0.00
## y1[3] 0.00
## y1[4] 1.00
## y1[5] 0.00
## y1[6] 3.00
## y1[7] 0.00
## y1[8] 0.00
## y1[9] 1.00
## y1[10] 1.00
## y1[11] 4.00
## y1[12] 2.00
## y1[13] 8.00
## y1[14] 0.00
## y1[15] 2.00
## y1[16] 2.00
## y1[17] 1.00
## y1[18] 0.00
## y1[19] 3.00
## y1[20] 0.00
## y1[21] 5.00
## y1[22] 8.00
## y1[23] 0.00
## y1[24] 0.00
## y1[25] 1.00
## y1[26] 1.00
## y1[27] 0.00
## y1[28] 0.00
## y1[29] 4.00
## y1[30] 1.00
## l1[1] 0.42
## l1[2] 0.42
## l1[3] -0.97
## l1[4] 0.10
## l1[5] 0.21
## l1[6] 0.62
## l1[7] 0.35
## l1[8] 0.38
## l1[9] -0.93
## l1[10] 0.72
## l1[11] 1.12
## l1[12] 1.16
## l1[13] 1.90
## l1[14] -1.08
## l1[15] 0.38
## l1[16] -0.64
## l1[17] -0.43
## l1[18] -0.93
## l1[19] 1.33
## l1[20] -0.75
## l1[21] 1.84
## l1[22] 1.39
## l1[23] -1.12
## l1[24] 0.21
## l1[25] -0.70
## l1[26] -0.33
## l1[27] -0.55
## l1[28] -0.33
## l1[29] 0.88
## l1[30] -0.32mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='l1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.44 -6.14 1.19 0.96 -8.80 -5.15 1.00 721 1058
## b[1] -0.26 -0.25 0.27 0.27 -0.72 0.14 1.00 916 905
## b[2] 1.65 1.64 0.29 0.30 1.18 2.14 1.00 1300 1192
## b[3] 0.55 0.54 0.29 0.28 0.08 1.03 1.00 888 839
## y1[1] 1.52 1.00 1.26 1.48 0.00 4.00 1.00 2041 1804
## y1[2] 1.54 1.00 1.28 1.48 0.00 4.00 1.00 1608 1747
## y1[3] 0.38 0.00 0.62 0.00 0.00 2.00 1.00 1975 1828
## y1[4] 1.10 1.00 1.12 1.48 0.00 3.00 1.00 2011 1763
## y1[5] 1.19 1.00 1.15 1.48 0.00 3.00 1.00 1831 1824
## y1[6] 1.81 2.00 1.37 1.48 0.00 4.00 1.00 1930 1896
## y1[7] 1.45 1.00 1.24 1.48 0.00 4.00 1.00 1677 1903
## y1[8] 1.45 1.00 1.22 1.48 0.00 4.00 1.00 2029 2014
## y1[9] 0.40 0.00 0.65 0.00 0.00 2.00 1.00 1897 1906
## y1[10] 2.03 2.00 1.47 1.48 0.00 5.00 1.00 2148 2108
## y1[11] 3.07 3.00 1.91 1.48 0.00 6.00 1.00 1960 1989
## y1[12] 3.17 3.00 1.84 1.48 1.00 6.00 1.00 1772 1858
## y1[13] 6.60 6.00 3.00 2.97 2.00 12.00 1.00 1955 1832
## y1[14] 0.34 0.00 0.59 0.00 0.00 1.00 1.00 1769 1797
## y1[15] 1.48 1.00 1.25 1.48 0.00 4.00 1.01 1866 1867
## y1[16] 0.52 0.00 0.74 0.00 0.00 2.00 1.00 1877 1938
## y1[17] 0.66 0.00 0.84 0.00 0.00 2.00 1.00 1860 1901
## y1[18] 0.42 0.00 0.67 0.00 0.00 2.00 1.00 1954 2006
## y1[19] 3.76 3.00 2.22 1.48 1.00 8.00 1.00 1874 1655
## y1[20] 0.47 0.00 0.70 0.00 0.00 2.00 1.00 2015 1920
## y1[21] 6.35 6.00 2.79 2.97 2.00 11.00 1.00 1780 1997
## y1[22] 3.93 4.00 2.18 1.48 1.00 8.00 1.00 2083 1962
## y1[23] 0.33 0.00 0.60 0.00 0.00 1.00 1.00 1948 2002
## y1[24] 1.22 1.00 1.12 1.48 0.00 3.00 1.00 1819 1720
## y1[25] 0.52 0.00 0.77 0.00 0.00 2.00 1.00 2107 1797
## y1[26] 0.72 0.00 0.89 0.00 0.00 2.00 1.00 1792 1714
## y1[27] 0.58 0.00 0.82 0.00 0.00 2.00 1.00 1991 1883
## y1[28] 0.70 1.00 0.84 1.48 0.00 2.00 1.00 1805 1648
## y1[29] 2.40 2.00 1.63 1.48 0.00 5.00 1.00 1770 1922
## y1[30] 0.74 1.00 0.87 1.48 0.00 2.00 1.00 1939 1882
## l1[1] 0.38 0.38 0.21 0.22 0.02 0.72 1.00 1591 1521
## l1[2] 0.39 0.39 0.23 0.23 -0.02 0.76 1.00 989 1077
## l1[3] -1.02 -1.00 0.35 0.36 -1.60 -0.48 1.00 949 851
## l1[4] 0.06 0.08 0.24 0.24 -0.35 0.44 1.00 935 984
## l1[5] 0.17 0.18 0.24 0.24 -0.23 0.54 1.00 951 948
## l1[6] 0.58 0.58 0.20 0.20 0.26 0.89 1.00 1698 1635
## l1[7] 0.32 0.33 0.23 0.23 -0.08 0.69 1.00 976 1058
## l1[8] 0.35 0.35 0.22 0.22 -0.02 0.69 1.00 1573 1546
## l1[9] -0.97 -0.96 0.35 0.36 -1.55 -0.44 1.00 946 820
## l1[10] 0.69 0.69 0.19 0.19 0.38 0.98 1.00 1758 1596
## l1[11] 1.09 1.09 0.26 0.25 0.64 1.50 1.00 1183 1266
## l1[12] 1.13 1.13 0.17 0.17 0.84 1.40 1.00 1862 1473
## l1[13] 1.87 1.88 0.22 0.22 1.52 2.22 1.00 1542 1619
## l1[14] -1.12 -1.10 0.37 0.38 -1.73 -0.56 1.00 959 874
## l1[15] 0.35 0.35 0.22 0.22 -0.03 0.69 1.00 1573 1546
## l1[16] -0.68 -0.67 0.31 0.32 -1.20 -0.20 1.00 922 806
## l1[17] -0.47 -0.46 0.29 0.29 -0.96 -0.03 1.00 908 798
## l1[18] -0.98 -0.97 0.41 0.43 -1.65 -0.33 1.01 1345 1202
## l1[19] 1.30 1.31 0.28 0.27 0.82 1.74 1.01 1204 1380
## l1[20] -0.79 -0.78 0.32 0.33 -1.34 -0.30 1.00 930 801
## l1[21] 1.81 1.82 0.21 0.21 1.46 2.15 1.00 1577 1628
## l1[22] 1.36 1.36 0.18 0.18 1.06 1.64 1.00 1816 1593
## l1[23] -1.17 -1.15 0.44 0.46 -1.89 -0.47 1.01 1337 1249
## l1[24] 0.18 0.18 0.24 0.24 -0.23 0.56 1.00 1510 1472
## l1[25] -0.74 -0.73 0.37 0.38 -1.35 -0.15 1.01 1362 1197
## l1[26] -0.37 -0.36 0.31 0.32 -0.90 0.14 1.01 1398 1367
## l1[27] -0.59 -0.58 0.35 0.36 -1.17 -0.03 1.01 1373 1289
## l1[28] -0.37 -0.37 0.28 0.28 -0.84 0.05 1.00 906 905
## l1[29] 0.84 0.85 0.24 0.24 0.42 1.23 1.00 1128 1179
## l1[30] -0.36 -0.35 0.31 0.32 -0.88 0.14 1.01 1400 1367
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3 4 6
## 0 7 3 0 0 0 0
## 1 4 3 0 0 0 0
## 2 0 3 1 0 0 0
## 3 0 1 1 1 0 1
## 4 0 0 1 1 1 0
## 5 0 0 0 1 0 0
## 7 0 0 0 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3 4 6
## 0 3 3 0 0 0 0
## 1 3 0 0 0 0 0
## 2 0 1 0 0 0 0
## 3 0 1 1 1 0 0
## 4 0 0 0 1 0 0
## 5 0 0 0 0 0 0
## 7 0 0 0 0 0 0
##
## , , = b
##
##
## 0 1 2 3 4 6
## 0 4 0 0 0 0 0
## 1 1 3 0 0 0 0
## 2 0 2 1 0 0 0
## 3 0 0 0 0 0 1
## 4 0 0 1 0 1 0
## 5 0 0 0 1 0 0
## 7 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 6
## 0 9 1 0 0 0
## 1 5 2 0 0 0
## 2 0 4 0 0 0
## 3 0 1 1 1 1
## 4 0 1 0 2 0
## 5 0 0 0 1 0
## 7 0 0 0 0 1cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 6
## 0 5 1 0 0 0
## 1 3 0 0 0 0
## 2 0 1 0 0 0
## 3 0 1 1 1 0
## 4 0 0 0 1 0
## 5 0 0 0 0 0
## 7 0 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3 6
## 0 4 0 0 0 0
## 1 2 2 0 0 0
## 2 0 3 0 0 0
## 3 0 0 0 0 1
## 4 0 1 0 1 0
## 5 0 0 0 1 0
## 7 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
logistic regression
# of samples is N,
objective variable 0/1 binary
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~Ber(pi), 0/1 binary
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-2.stan
logistic regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~bernoulli_logit(z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=bernoulli_rng(inv_logit(z1[i]));
}
}n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='binomial') # it can caluculte when all trials are once##
## Call: glm(formula = f, family = "binomial", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.201 0.167 0.324
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 41.1
## Residual Deviance: 40.9 AIC: 46.9X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -38.9076
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -20.4263 0.000312115 7.48949e-05 0.8677 0.8677 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -20.43
## b[1] 0.20
## b[2] 0.17
## b[3] 0.32
## y1[1] 1.00
## y1[2] 0.00
## y1[3] 1.00
## y1[4] 1.00
## y1[5] 0.00
## y1[6] 0.00
## y1[7] 1.00
## y1[8] 1.00
## y1[9] 1.00
## y1[10] 1.00
## y1[11] 1.00
## y1[12] 0.00
## y1[13] 1.00
## y1[14] 0.00
## y1[15] 0.00
## y1[16] 1.00
## y1[17] 1.00
## y1[18] 1.00
## y1[19] 0.00
## y1[20] 1.00
## y1[21] 1.00
## y1[22] 1.00
## y1[23] 1.00
## y1[24] 1.00
## y1[25] 1.00
## y1[26] 1.00
## y1[27] 0.00
## y1[28] 1.00
## y1[29] 1.00
## y1[30] 0.00
## z1[1] 0.35
## z1[2] 0.18
## z1[3] 0.53
## z1[4] 0.17
## z1[5] 0.10
## z1[6] 0.27
## z1[7] 0.16
## z1[8] 0.09
## z1[9] 0.31
## z1[10] 0.21
## z1[11] 0.07
## z1[12] 0.23
## z1[13] 0.59
## z1[14] 0.62
## z1[15] 0.12
## z1[16] 0.52
## z1[17] 0.28
## z1[18] 0.12
## z1[19] 0.19
## z1[20] 0.20
## z1[21] 0.42
## z1[22] 0.56
## z1[23] 0.15
## z1[24] 0.40
## z1[25] 0.08
## z1[26] 0.15
## z1[27] 0.13
## z1[28] 0.11
## z1[29] 0.37
## z1[30] 0.46mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.00 -21.68 1.26 1.03 -24.55 -20.61 1.00 917 1112
## b[1] 0.22 0.23 0.47 0.46 -0.54 1.02 1.00 1639 1365
## b[2] 0.17 0.13 0.81 0.78 -1.13 1.56 1.00 1388 1220
## b[3] 0.36 0.36 0.91 0.91 -1.07 1.88 1.00 1083 970
## y1[1] 0.58 1.00 0.49 0.00 0.00 1.00 1.00 1846 NA
## y1[2] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 2016 NA
## y1[3] 0.64 1.00 0.48 0.00 0.00 1.00 1.00 1982 NA
## y1[4] 0.54 1.00 0.50 0.00 0.00 1.00 1.00 1859 NA
## y1[5] 0.52 1.00 0.50 0.00 0.00 1.00 1.00 1963 NA
## y1[6] 0.57 1.00 0.50 0.00 0.00 1.00 1.00 2047 NA
## y1[7] 0.55 1.00 0.50 0.00 0.00 1.00 1.00 1954 NA
## y1[8] 0.52 1.00 0.50 0.00 0.00 1.00 1.00 1958 NA
## y1[9] 0.56 1.00 0.50 0.00 0.00 1.00 1.00 1820 NA
## y1[10] 0.56 1.00 0.50 0.00 0.00 1.00 1.00 1874 NA
## y1[11] 0.51 1.00 0.50 0.00 0.00 1.00 1.00 1939 NA
## y1[12] 0.57 1.00 0.50 0.00 0.00 1.00 1.00 1964 NA
## y1[13] 0.64 1.00 0.48 0.00 0.00 1.00 1.00 1811 NA
## y1[14] 0.65 1.00 0.48 0.00 0.00 1.00 1.00 1733 NA
## y1[15] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 1855 NA
## y1[16] 0.62 1.00 0.48 0.00 0.00 1.00 1.00 1833 NA
## y1[17] 0.55 1.00 0.50 0.00 0.00 1.00 1.00 1841 NA
## y1[18] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 1973 NA
## y1[19] 0.55 1.00 0.50 0.00 0.00 1.00 1.00 1932 NA
## y1[20] 0.56 1.00 0.50 0.00 0.00 1.00 1.00 1894 NA
## y1[21] 0.60 1.00 0.49 0.00 0.00 1.00 1.00 1879 NA
## y1[22] 0.63 1.00 0.48 0.00 0.00 1.00 1.00 1746 NA
## y1[23] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 1934 NA
## y1[24] 0.61 1.00 0.49 0.00 0.00 1.00 1.00 2047 NA
## y1[25] 0.52 1.00 0.50 0.00 0.00 1.00 1.00 2039 NA
## y1[26] 0.52 1.00 0.50 0.00 0.00 1.00 1.00 1605 NA
## y1[27] 0.56 1.00 0.50 0.00 0.00 1.00 1.00 2089 NA
## y1[28] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 2022 NA
## y1[29] 0.59 1.00 0.49 0.00 0.00 1.00 1.00 1847 NA
## y1[30] 0.62 1.00 0.48 0.00 0.00 1.00 1.00 1954 NA
## z1[1] 0.37 0.35 0.95 0.94 -1.12 1.97 1.00 1389 1219
## z1[2] 0.20 0.20 0.45 0.43 -0.55 0.95 1.00 1711 1389
## z1[3] 0.59 0.57 0.77 0.75 -0.64 1.90 1.00 1401 1542
## z1[4] 0.18 0.19 0.45 0.43 -0.57 0.93 1.00 1731 1330
## z1[5] 0.11 0.13 0.58 0.58 -0.81 1.08 1.00 1670 1518
## z1[6] 0.28 0.28 0.63 0.61 -0.70 1.31 1.00 1471 1328
## z1[7] 0.17 0.18 0.46 0.44 -0.58 0.94 1.00 1747 1366
## z1[8] 0.10 0.12 0.62 0.62 -0.87 1.12 1.00 1636 1518
## z1[9] 0.33 0.30 0.79 0.77 -0.90 1.64 1.00 1424 1340
## z1[10] 0.23 0.25 0.49 0.48 -0.57 1.06 1.00 1596 1519
## z1[11] 0.08 0.10 0.68 0.69 -1.00 1.21 1.00 1580 1536
## z1[12] 0.25 0.26 0.52 0.52 -0.60 1.13 1.00 1550 1406
## z1[13] 0.65 0.63 0.85 0.82 -0.72 2.05 1.00 1403 1370
## z1[14] 0.68 0.68 0.92 0.92 -0.78 2.14 1.00 1399 1364
## z1[15] 0.13 0.14 0.53 0.52 -0.71 1.00 1.00 1709 1471
## z1[16] 0.58 0.56 0.77 0.75 -0.63 1.89 1.00 1402 1335
## z1[17] 0.30 0.28 0.68 0.66 -0.77 1.42 1.00 1449 1320
## z1[18] 0.13 0.14 0.53 0.52 -0.71 1.00 1.00 1709 1471
## z1[19] 0.21 0.22 0.46 0.44 -0.55 0.98 1.00 1672 1430
## z1[20] 0.22 0.23 0.47 0.46 -0.56 1.03 1.00 1628 1449
## z1[21] 0.47 0.43 0.89 0.89 -0.98 1.95 1.00 1410 1365
## z1[22] 0.62 0.59 0.81 0.79 -0.65 1.95 1.00 1401 1501
## z1[23] 0.17 0.17 0.47 0.45 -0.59 0.93 1.00 1751 1406
## z1[24] 0.45 0.42 0.95 0.95 -1.08 2.01 1.00 1409 1322
## z1[25] 0.09 0.11 0.66 0.66 -0.95 1.18 1.00 1600 1540
## z1[26] 0.16 0.17 0.47 0.45 -0.60 0.94 1.00 1753 1406
## z1[27] 0.14 0.14 0.51 0.49 -0.68 0.97 1.00 1723 1519
## z1[28] 0.13 0.13 0.55 0.54 -0.76 1.04 1.00 1692 1448
## z1[29] 0.39 0.35 1.02 1.01 -1.24 2.11 1.00 1383 1203
## z1[30] 0.52 0.48 0.80 0.79 -0.78 1.85 1.00 1402 1405
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 1
## 0 13
## 1 17cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 1
## 0 10
## 1 12
##
## , , = b
##
##
## 1
## 0 3
## 1 5par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 1
## 0 13
## 1 17cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 1
## 0 10
## 1 12
##
## , , = b
##
## map
## 1
## 0 3
## 1 5par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
multi logistic regression
# of samples is N,
# of trials of a sample i is mi,
objective variable [0,n], integer
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~B(mi, pi), # of acutual incident
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-3.stan
multi logistic regression
data{
int N;
int K;
array[N] int m;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~binomial_logit(m,z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
}
}n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)
mdl=cmdstan_model('./ex6-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -155.274
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -96.8953 0.000390871 8.9511e-05 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -96.90
## b[1] 0.07
## b[2] 0.88
## b[3] 0.99
## y1[1] 2.00
## y1[2] 2.00
## y1[3] 0.00
## y1[4] 3.00
## y1[5] 3.00
## y1[6] 0.00
## y1[7] 7.00
## y1[8] 4.00
## y1[9] 5.00
## y1[10] 5.00
## y1[11] 1.00
## y1[12] 6.00
## y1[13] 7.00
## y1[14] 2.00
## y1[15] 2.00
## y1[16] 4.00
## y1[17] 0.00
## y1[18] 2.00
## y1[19] 1.00
## y1[20] 2.00
## y1[21] 4.00
## y1[22] 3.00
## y1[23] 6.00
## y1[24] 4.00
## y1[25] 3.00
## y1[26] 2.00
## y1[27] 3.00
## y1[28] 3.00
## y1[29] 4.00
## y1[30] 2.00
## z1[1] 1.48
## z1[2] 0.84
## z1[3] -0.34
## z1[4] 1.83
## z1[5] -0.34
## z1[6] -0.50
## z1[7] 1.38
## z1[8] 0.51
## z1[9] 1.29
## z1[10] -0.58
## z1[11] -0.40
## z1[12] 0.65
## z1[13] 1.87
## z1[14] -0.67
## z1[15] -0.16
## z1[16] 0.42
## z1[17] -0.71
## z1[18] 0.31
## z1[19] 0.29
## z1[20] 0.92
## z1[21] 0.29
## z1[22] 0.51
## z1[23] 1.81
## z1[24] 0.90
## z1[25] 1.22
## z1[26] -0.73
## z1[27] -0.67
## z1[28] 1.03
## z1[29] 0.58
## z1[30] 0.53mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -98.46 -98.11 1.32 1.03 -101.16 -97.07 1.00 694 1140
## b[1] 0.08 0.07 0.33 0.33 -0.45 0.65 1.00 624 925
## b[2] 0.91 0.91 0.37 0.36 0.32 1.54 1.00 908 1081
## b[3] 1.01 1.01 0.37 0.37 0.42 1.62 1.00 682 943
## y1[1] 2.43 3.00 0.69 0.00 1.00 3.00 1.00 1994 NA
## y1[2] 1.39 1.00 0.66 1.48 0.00 2.00 1.00 2110 NA
## y1[3] 0.41 0.00 0.49 0.00 0.00 1.00 1.00 1849 NA
## y1[4] 2.57 3.00 0.62 0.00 1.00 3.00 1.00 1787 NA
## y1[5] 2.51 2.00 1.25 1.48 1.00 5.00 1.00 1893 1772
## y1[6] 0.74 1.00 0.69 1.48 0.00 2.00 1.00 1912 NA
## y1[7] 7.19 7.00 1.27 1.48 5.00 9.00 1.00 1844 NA
## y1[8] 4.36 4.00 1.40 1.48 2.00 7.00 1.00 1520 NA
## y1[9] 4.73 5.00 1.02 1.48 3.00 6.00 1.00 1924 NA
## y1[10] 3.22 3.00 1.57 1.48 1.00 6.00 1.00 1646 1746
## y1[11] 1.17 1.00 0.87 1.48 0.00 3.00 1.00 2021 NA
## y1[12] 4.64 5.00 1.28 1.48 3.00 7.00 1.00 1970 NA
## y1[13] 6.89 7.00 1.07 1.48 5.00 8.00 1.00 1882 NA
## y1[14] 3.06 3.00 1.52 1.48 1.00 6.00 1.00 1823 1814
## y1[15] 4.18 4.00 1.67 1.48 2.00 7.00 1.00 1479 1832
## y1[16] 5.49 6.00 1.58 1.48 3.00 8.00 1.00 1942 2088
## y1[17] 0.31 0.00 0.46 0.00 0.00 1.00 1.00 1962 NA
## y1[18] 1.14 1.00 0.71 1.48 0.00 2.00 1.00 1893 NA
## y1[19] 1.72 2.00 0.87 1.48 0.00 3.00 1.00 1953 NA
## y1[20] 1.43 2.00 0.65 0.00 0.00 2.00 1.00 2025 NA
## y1[21] 2.80 3.00 1.18 1.48 1.00 5.00 1.00 1453 NA
## y1[22] 3.75 4.00 1.23 1.48 2.00 6.00 1.00 1765 NA
## y1[23] 5.17 5.00 0.88 1.48 4.00 6.00 1.00 1708 NA
## y1[24] 4.27 4.00 1.11 1.48 2.00 6.00 1.00 1816 NA
## y1[25] 2.30 2.00 0.75 1.48 1.00 3.00 1.00 1876 NA
## y1[26] 2.00 2.00 1.22 1.48 0.00 4.00 1.00 1931 1848
## y1[27] 1.98 2.00 1.21 1.48 0.00 4.00 1.00 1947 1932
## y1[28] 5.20 5.00 1.17 1.48 3.00 7.00 1.00 1846 NA
## y1[29] 5.16 5.00 1.44 1.48 3.00 7.00 1.00 2117 1799
## y1[30] 3.80 4.00 1.23 1.48 2.00 6.00 1.00 2144 NA
## z1[1] 1.52 1.52 0.36 0.34 0.95 2.10 1.00 1238 1102
## z1[2] 0.86 0.85 0.23 0.22 0.50 1.25 1.00 1704 1540
## z1[3] -0.35 -0.35 0.28 0.27 -0.79 0.12 1.00 780 821
## z1[4] 1.89 1.88 0.48 0.46 1.11 2.68 1.00 1110 1151
## z1[5] -0.34 -0.35 0.28 0.27 -0.78 0.13 1.00 779 821
## z1[6] -0.51 -0.51 0.28 0.27 -0.96 -0.02 1.00 912 982
## z1[7] 1.42 1.42 0.33 0.31 0.89 1.95 1.00 1295 1185
## z1[8] 0.53 0.51 0.46 0.46 -0.21 1.33 1.00 626 809
## z1[9] 1.33 1.33 0.30 0.29 0.84 1.82 1.00 1358 1248
## z1[10] -0.59 -0.59 0.29 0.28 -1.06 -0.11 1.00 982 1190
## z1[11] -0.40 -0.41 0.28 0.27 -0.85 0.07 1.00 823 947
## z1[12] 0.66 0.65 0.24 0.23 0.28 1.06 1.00 1465 1598
## z1[13] 1.92 1.92 0.49 0.47 1.12 2.74 1.00 1101 1151
## z1[14] -0.68 -0.68 0.30 0.30 -1.18 -0.19 1.00 1067 1199
## z1[15] -0.16 -0.16 0.29 0.29 -0.62 0.34 1.00 682 747
## z1[16] 0.42 0.41 0.28 0.28 -0.03 0.89 1.00 1180 1373
## z1[17] -0.73 -0.73 0.31 0.30 -1.24 -0.21 1.00 1104 1153
## z1[18] 0.32 0.31 0.31 0.30 -0.19 0.84 1.00 1092 1485
## z1[19] 0.29 0.28 0.32 0.31 -0.23 0.81 1.00 1074 1374
## z1[20] 0.94 0.93 0.23 0.22 0.57 1.32 1.00 1698 1426
## z1[21] 0.29 0.28 0.31 0.31 -0.22 0.82 1.00 1078 1439
## z1[22] 0.52 0.51 0.26 0.26 0.10 0.96 1.00 1284 1654
## z1[23] 1.87 1.86 0.47 0.45 1.11 2.65 1.00 1117 1165
## z1[24] 0.92 0.92 0.23 0.22 0.55 1.31 1.00 1702 1450
## z1[25] 1.25 1.25 0.28 0.27 0.79 1.71 1.00 1425 1254
## z1[26] -0.75 -0.75 0.31 0.31 -1.27 -0.22 1.00 1117 1111
## z1[27] -0.69 -0.69 0.30 0.30 -1.19 -0.19 1.00 1074 1198
## z1[28] 1.06 1.06 0.24 0.23 0.67 1.47 1.00 1600 1223
## z1[29] 0.59 0.58 0.25 0.24 0.19 1.01 1.00 1371 1599
## z1[30] 0.54 0.52 0.26 0.25 0.13 0.97 1.00 1305 1626
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3 4 5 6 7
## 0 1 2 0 0 0 0 0 0
## 1 1 2 2 0 0 0 0 0
## 2 0 0 3 2 0 0 0 0
## 3 0 0 0 2 0 1 0 0
## 4 0 0 1 1 3 0 0 0
## 5 0 0 0 0 1 4 0 0
## 6 0 0 0 0 1 0 1 0
## 7 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3 4 5 6 7
## 0 1 0 0 0 0 0 0 0
## 1 1 2 0 0 0 0 0 0
## 2 0 0 2 2 0 0 0 0
## 3 0 0 0 0 0 0 0 0
## 4 0 0 1 0 1 0 0 0
## 5 0 0 0 0 1 0 0 0
## 6 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0
##
## , , = b
##
##
## 0 1 2 3 4 5 6 7
## 0 0 2 0 0 0 0 0 0
## 1 0 0 2 0 0 0 0 0
## 2 0 0 1 0 0 0 0 0
## 3 0 0 0 2 0 1 0 0
## 4 0 0 0 1 2 0 0 0
## 5 0 0 0 0 0 4 0 0
## 6 0 0 0 0 1 0 1 0
## 7 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 4 5 6 7 8
## 0 1 1 1 0 0 0 0 0 0
## 1 1 2 2 0 0 0 0 0 0
## 2 0 0 1 4 0 0 0 0 0
## 3 0 0 0 2 0 0 1 0 0
## 4 0 0 1 1 1 2 0 0 0
## 5 0 0 0 0 1 2 2 0 0
## 6 0 0 0 0 1 0 1 0 0
## 7 0 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 1 0cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 4 5 6 7 8
## 0 1 0 0 0 0 0 0 0 0
## 1 1 2 0 0 0 0 0 0 0
## 2 0 0 1 3 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0
## 4 0 0 1 0 0 1 0 0 0
## 5 0 0 0 0 1 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3 4 5 6 7 8
## 0 0 1 1 0 0 0 0 0 0
## 1 0 0 2 0 0 0 0 0 0
## 2 0 0 0 1 0 0 0 0 0
## 3 0 0 0 2 0 0 1 0 0
## 4 0 0 0 1 1 1 0 0 0
## 5 0 0 0 0 0 2 2 0 0
## 6 0 0 0 0 1 0 1 0 0
## 7 0 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 1 0par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
gamma regression
objective variable [0,Infinity)
# of samples is N,
log (a/li)=sum(bj*xji),j=[0,K],i=[1,N]
li=a/exp(sum(bj*xji))
yi~Ga(a,li) ex6-4.stan
gamma regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> a;
}
model{
vector[N] l;
for(i in 1:N){
l[i]=a/exp(X[i]*b);
}
y~gamma(a,l);
}n=20
tb=tibble(x1=runif(n,0,2),x2=runif(n,0,2),
y=rgamma(n,3,3/exp(x1+x2)))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-4.stan')
mle=mdl$optimize(data=data) # it sometimes occur error and stop process## Initial log joint probability = -839.095
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -48.9895 0.000130897 0.000166455 1 1 20
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -48.99
## b[1] -0.09
## b[2] 1.27
## b[3] 0.88
## a 3.59mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -49.92 -49.57 1.59 1.34 -53.02 -48.04 1.00 770 1126
## b[1] -0.07 -0.08 0.30 0.30 -0.58 0.41 1.00 1002 1010
## b[2] 1.28 1.28 0.27 0.25 0.83 1.73 1.00 1082 1105
## b[3] 0.88 0.87 0.25 0.25 0.49 1.29 1.00 1315 1099
## a 3.46 3.34 1.12 1.03 1.88 5.50 1.00 932 876
negative binomial regression
The event with probability p do not occur y times till the event occur a times
(negative binomial distribution has larger variance than poisson distribution)
y~NB(a,p), log(p)=X*b
y~NB(a,l0/(1+l0)) when y~Po(l), l~Ga(a,l0), l0=a/E[l]ex6-5.stan
negative binomial regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> a;
}
model{
a~cauchy(0,5);
y~neg_binomial_2_log(X*b,a);
}n=20
tb=tibble(x1=runif(n,-1,0),x2=runif(n,-1,0),
y=rnbinom(n,3,exp(x1+x2)))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-5.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -66.829
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -46.4879 0.0203844 0.000489725 1 1 19
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -46.49
## b[1] -0.19
## b[2] -1.28
## b[3] -1.63
## a 3.20mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -47.32 -46.94 1.52 1.28 -50.26 -45.59 1.01 647 710
## b[1] -0.20 -0.16 0.53 0.52 -1.11 0.60 1.01 743 576
## b[2] -1.30 -1.29 0.62 0.60 -2.34 -0.32 1.01 1023 1096
## b[3] -1.62 -1.60 0.76 0.76 -2.87 -0.44 1.00 910 798
## a 4.37 3.49 3.40 1.79 1.52 9.38 1.00 1134 826
beta regression
using prior of binomial distribution parameter p[0,1]
y~B(n,p), p~beta(a,b)
m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)
m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*sex6-6.stan
beta regression
data{
int N;
int K;
vector[N] p;
matrix[N,K] X;
}
parameters{
vector[K] be;
real<lower=0> s;
}
model{
vector[N] a;
vector[N] b;
for(i in 1:N){
a[i]=inv_logit(X[i]*be)*s;
b[i]=(1-inv_logit(X[i]*be))*s;
}
p~beta(a,b);
}n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$p)
plot(tb$x2,tb$p)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),p=tb$p,X=X)
mdl=cmdstan_model('./ex6-6.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -23.3907
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 5.89803 0.000126732 0.00013279 1 1 28
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 5.90
## be[1] -1.98
## be[2] 4.66
## be[3] 1.80
## s 4.56mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 5.37 5.73 1.48 1.21 2.38 7.09 1.01 610 867
## be[1] -2.07 -2.06 0.68 0.68 -3.18 -0.93 1.01 636 804
## be[2] 4.89 4.86 1.56 1.57 2.34 7.43 1.01 726 631
## be[3] 1.84 1.87 1.73 1.73 -0.95 4.59 1.00 656 643
## s 4.64 4.50 1.37 1.30 2.65 7.12 1.00 1644 1434
beta binomial regression
fitting to distribution has larger variance than binomial distribution
y~betaB(n,a,b) when y~B(n,p), p~beta(a,b)
m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)
m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*sex6-7.stan
beta binomial regression
data{
int N;
int K;
array[N] int y;
array[N] int n;
matrix[N,K] X;
}
parameters{
vector[K] be;
real<lower=0> s;
}
model{
vector[N] a;
vector[N] b;
for(i in 1:N){
a[i]=inv_logit(X[i]*be)*s;
b[i]=(1-inv_logit(X[i]*be))*s;
}
y~beta_binomial(n,a,b);
s~cauchy(0,5);
}n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3),
n1=floor(runif(n,5,9)),
y=rbinom(n,n1,p))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),n=tb$n1,y=tb$y,X=X)
mdl=cmdstan_model('./ex6-7.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -80.2502
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -78.1212 0.000434371 0.00119507 0.5333 0.5333 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -78.12
## be[1] -1.73
## be[2] 3.87
## be[3] 2.41
## s 3.38mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -78.84 -78.56 1.41 1.31 -81.61 -77.14 1.00 735 1078
## be[1] -1.95 -1.92 0.78 0.73 -3.33 -0.75 1.00 919 1085
## be[2] 4.17 4.20 1.88 1.89 1.05 7.27 1.00 1003 958
## be[3] 2.86 2.76 2.17 2.13 -0.59 6.64 1.00 853 866
## s 5.38 4.06 9.63 2.06 1.81 11.55 1.01 1944 1042
